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Nonequilibrium  phenomena  are  well  understood  for  states  near  equilibrium. 
However,  even  the  qualitative  properties  of  transport  and  fluctuations  far  from 
equilibrium  are  poorly  understood  due  to  the  complexity  of  the  physical  states. 
The  objective  here  is  to  explore  in  detail  dynamical  properties  of  a  simple  fluid 
far  from  equilibrium  for  the  specific  case  of  uniform  shear  flow.  The  analysis  is 
limited  to  conditions  of  low  density  for  which  there  is  a  controlled  kinetic  the- 
ory description  based  in  the  principles  of  nonequilibrium  statistical  mechanics. 
These  kinetic  equations  provide  a  description  of  both  transport  and  fluctuations 
in  general  nonequilibrium  states.  However  it  is  very  difficult  to  solve  these  equa- 
tions. To  make  this  description  practical  far  from  equilibrium,  a  related  set  of 
model  kinetic  equations  is  constructed  by  retaining  the  essential  mathematical 
and  physical  properties  of  the  original  kinetic  equations,  while  admitting  more 


practical  analysis.  These  model  equations  are  solved  for  the  reference  state  of 
uniform  shear  flow,  fluctuations  about  that  state,  and  for  the  dynamics  of  small 
perturbations  from  that  state.  In  this  shear  flow,  the  shear  rate  is  the  only  con- 
trol parameter  that  measures  how  far  the  system  is  from  equilibrium.  First,  the 
distribution  function  for  the  reference  state  is  calculated  exactly  and  an  expan- 
sion in  small  spatial  deviations  from  the  reference  state  is  performed  to  describe 
the  dynamics  of  states  near  uniform  shear  flow.  This  solution  to  the  kinetic 
equation  also  provides  a  hydrodynamic  description,  and  a  closed  set  of  gener- 
alized hydrodynamic  equations  is  derived.  The  associated  nonlinear  transport 
coefficients  are  identified  and  calculated  for  an  arbitrary  shear  rate.  The  hydro- 
dynamic  modes  defined  by  the  linearized  hydrodynamic  equations  are  studied  to 
determine  the  shear  rate  dependence  of  relaxation  back  to  the  reference  state. 
A  new  long  wavelength  instability  for  uniform  shear  flow  is  found  such  that  the 
hydrodynamic  modes  are  growing  initially  for  finite  shear  rate  and  sufficiently 
long  wavelength.  Finally,  spatially  long  ranged  correlations  within  the  uniform 
shear  state  are  studied  at  large  shear  rates. 


CHAPTER  1 
INTRODUCTION 

The  world  in  which  we  live  consists  of  a  large  number  of  particles  which  inter- 
act in  a  very  complex  way.  Statistical  mechanics  is  the  field  of  describing  static 
and  dynamical  processes  of  many  body  systems,  using  probabilistic  concepts. 
The  state  is  specified  not  by  the  position  and  momenta  of  individual  particles  (it 
is  not  wrong,  only  impractical),  but  by  a  distribution  function,  p(t),  which  is  a 
probability  density  of  the  system  in  the  N-particle  phase  space.  The  dynamics 
of  this  distribution  function  is  calculated  from  the  Liouville  equation.  With  the 
knowledge  of  this  function,  the  average  value  of  any  dynamical  variable  (physical 
observable)  can  be  obtained  in  principle.  However,  these  average  values  are  not 
meaningful  for  describing  the  system  unless  the  fluctuations  about  the  average 
and  actual  values  are  known.  These  fluctuations  also  are  related  closely  to  a 
number  of  experimental  methods  used  to  probe  many  body  systems  (e.g.  elastic 
and  inelastic  scattering  of  particles  and  radiation).  Therefore,  fluctuations  are 
very  important  properties  to  understand.  The  most  general  type  of  average  and 
the  fluctuations  are  defined  by 

{A(t))t0=JdTA(r,t)p(T,t0),  (1.1) 

(8A(t)5B(t'))t0  =  J  dr5A{T,t)SB(r,f)p{V,t0)t  (1.2) 


where  A(t),  B(t')  correspond  to  physical  observables  of  interest  at  times  t  and 
£',  and  the  bracket  means  the  average  over  the  distribution  function  at  time  io 
when  the  system  is  prepared.  Also,  T  denotes  a  point  in  the  N-particle  phase 
space,  representing  the  positions  and  momenta  of  all  particles,  and  SA(t)  means 
A(t)  -  (A(t))t0. 

At  equilibrium,  the  properties  of  the  average  values  and  their  fluctuations 
are  well  understood,  and  extensive  methods  exist  for  their  calculations  and  mea- 
surement. For  instance,  the  averages  provide  the  thermodynamic  variables  (e.g. 
internal  energy),  while  the  susceptibilities  are  related  to  fluctuations  at  equal 
times  (e.g.  specific  heat).  Both  of  these  can  be  calculated  as  averages  over  the 
known  Gibbs  distribution  function  using  Eqs.(l.l)  and  (1.2).  The  dynamics  of 
fluctuations  at  different  times  in  the  equilibrium  state  also  is  important  since  it 
describes  the  linear  response  of  an  equilibrium  system  to  small  initial  or  external 
perturbations  [1].  These  fluctuations  also  are  well  studied  both  theoretically  and 
experimentally. 

In  nonequilibrium  states,  the  methods  for  obtaining  the  values  of  the  physical 
observables  and  their  fluctuations  (both  static  and  dynamic)  are  still  in  the  pro- 
cess of  developing.  Significant  progress  has  been  made  for  states  near  equilibrium. 
For  example,  kinetic  and  hydrodynamic  equations  are  accurate  and  practical  for 
calculating  average  properties  near  equilibrium.  Furthermore,  the  dynamics  of 
fluctuations  in  states  near  equilibrium  has  been  studied  recently  in  some  de- 
tail [2j.  These  studies  show  interesting  qualitative  differences  from  the  corre- 
sponding fluctuations  at  equilibrium,  such  as  long  ranged  spatial  correlations  [3]. 


Recent  accurate  light  scattering  experiments  [4]  support  these  theoretical  predic- 
tions. Also,  these  studies  have  provided  a  simple  phenomenological  description, 
e.g.  fluctuating  hydrodynamics  [5],  to  understand  the  mechanisms  responsible 
for  these  correlations.  However,  for  the  states  far  from  equilibrium,  technical  dif- 
ficulties prohibit  having  any  significant  qualitative  and  quantitative  knowledge  of 
either  average  values  or  their  fluctuations  in  general.  The  purpose  of  the  work  in 
this  thesis  is  to  choose  a  special  state  far  from  equilibrium  and  to  explore  all  pos- 
sible aspects  of  static  and  dynamic  fluctuations  presently  unknown.  Progress  on 
this  difficult  problem  has  been  possible  because  there  is  a  well-developed  kinetic 
theory  at  low  density  and  in  principle  those  kinetic  equations  are  solvable  for  any 
nonequilibrium  state.  However  their  application  is  limited  for  practical  purposes. 
In  this  thesis,  "model"  kinetic  equations  are  developed  for  static  and  dynamic 
correlation  functions  by  using  fundamental  functional  relationships  between  the 
related  kinetic  equations  for  transport  and  fluctuations.  These  equations  provide 
a  practical  basis  for  describing  phenomena  in  general  nonequilibrium  states.  To 
illustrate  the  practicality  of  these  kinetic  equations,  we  choose  uniform  shear  flow 
state  as  a  prototype  nonequilibrium  test  case.  The  uniform  shear  flow  state  is 
characterized  by  a  very  simple  macroscopic  state.  This  flow  is  defined  by  a  lin- 
ear velocity  profile  along  the  x-axis  with  a  constant  velocity  gradient  along  the 
y-axis,  and  constant  density  and  temperature, 

ul  =  aijrJ,   a.ij  =  a5iX8jy,  n  =  constant,   T  =  T(t).  (1.3) 

where  u,  n,  T  are  flow  velocity,  density,  temperature  respectively.   The  velocity 


gradient  or  shear  rate,  a,  is  constant  and  represents  the  single  control  parameter 
measuring  how  far  the  system  is  from  equilibrium.  If  we  perform  local  Galilean 
transformations  r'  =  r  -  u(r)t,  v'  =  v  -  u(r),  the  macroscopic  state  can  be 
transformed  into  the  uniform  local  rest  state.  Because  of  this  property,  uniform 
shear  flow  can  be  generated  easily  with  Lees-Edwards  boundary  conditions  [6] 
which  are  simple  periodic  boundary  conditions  in  the  local  rest  frame.  Also,  the 
shear  effects  in  uniform  shear  flow  induce  viscous  heating  and  cause  the  temper- 
ature to  increase  in  time.  This  heating  problem  can  be  avoided  by  introducing 
a  thermostat  such  that  the  state  becomes  stationary.  Because  of  this  stationar- 
ity,  simplicity  of  the  macroscopic  state,  and  the  ideal  boundary  conditions,  this 
state  has  been  used  as  a  prototype  for  the  studies  of  properties  in  nonequilibrium 
systems.  There  have  been  many  computer  simulations  [6-11]  for  transport  prop- 
erties for  this  state.  However,  theoretical  studies  have  been  restricted  primarily 
to  the  limit  of  small  shear  rate  [2,12-17].  In  this  thesis,  we  explore  properties 
such  as  nonlinear  transport  and  fluctuations  for  states  far  from  equilibrium  for 
uniform  shear  flow  and  states  near  the  uniform  shear  flow. 

Traditionally  there  are  two  ways  to  calculate  transport  properties  and  fluctua- 
tions in  nonequilibrium  states.  The  first  is  a  direct  statistical  mechanics  approach. 
In  this  method,  a  formally  exact  solution  to  the  Liouville  equation  is  constructed. 
Use  of  this  solution  in  (1.1)  and  (1.2)  gives  an  explicit  representation  of  the  phys- 
ical observables  and  their  fluctuations.  Generally,  however,  it  is  very  difficult  to 
obtain  this  solution  in  a  practical  form.  Instead,  a  perturbation  method  can  be 
used  if  the  state  is  not  far  from  equilibrium  to  construct  an  explicit  form  for  p 


as  a  small  deviation  from  the  Gibbs  distribution.  In  Chapter  2,  the  correlation 
function  near  the  equilibrium  state  is  represented  in  this  way  for  uniform  shear 
flow.  The  resulting  form  is  suitable  for  study  via  standard  equilibrium  molecular 
simulation,  but  still  difficult  to  calculate  theoretically.  However  at  low  density,  it 
is  shown  that  the  calculation  is  tractable  and  the  nonequilibrium  density  fluctu- 
ations are  calculated  as  an  illustration.  The  result  shows  the  spatially  long  range 
correlations  mentioned  above,  in  agreement  with  previous  phenomenological  cal- 
culations [18].  This  brief  analysis  is  presented  to  demonstrate  the  new  effects 
that  occur  for  nonequilibrium  fluctuations  and  to  show  the  practical  limitations 
of  the  direct  computational  approach  for  states  far  from  equilibrium. 

The  second  method  to  describe  transport  and  fluctuations  is  by  the  kinetic 
theory.  It  is  this  method  that  is  used  in  the  remainder  of  the  thesis.  The  average  of 
any  physical  observable  that  is  represented  by  a  sum  of  single  particle  functions 
may  be  obtained  as  an  integral  over  the  solution  to  the  single  particle  kinetic 
equation,  /"', 

{A(t))t0=  fdxfMfatoMx.t);   A(t)  =  Y,^,(x,t).  (1.4) 

:=1 

Traditional  kinetic  theory  provides  an  equation  for  the  single  particle  reduced 
distribution  function  derived  from  an  approximation  to  the  BBGKY  hierarchy. 
At  low  density,  the  relevant  approximation  leads  to  the  famous  Boltzmann  equa- 
tion. Closely  related  kinetic  equations  for  fluctuations  also  can  be  constructed 
by  applying  the  same  many  body  analysis  used  to  derive  the  Boltzmann  equa- 
tion [19].  Generally,  these  kinetic  equations  are  very  difficult  to  solve  except  for 


a  few  special  cases,  or  for  states  near  equilibrium.  Sophisticated  Monte  Carlo 
simulation  methods  are  available  to  solve  the  Boltzmann  equation  [20],  but  they 
have  not  been  adapted  to  solve  the  equations  for  fluctuations.  Exact  solutions  to 
the  Boltzmann  equation  far  from  equilibrium  are  still  very  few.  A  formal  solution 
exists  for  uniform  shear  flow  [21,22],  but  this  solution  provides  only  equations 
for  moments  of  the  distribution  function  rather  than  for  the  distribution  func- 
tion itself.  In  the  case  of  the  Boltzmann  equation,  these  difficulties  have  been 
circumvented  by  replacing  the  collision  operator  with  a  simpler  form.  The  result 
is  a  "model"  kinetic  equation  that  preserves  the  most  important  mathematical 
and  physical  properties  of  the  Boltzmann  equation,  while  allowing  more  prac- 
tical analysis.  The  choice  considered  here  is  the  well  known  BGK  model  [23]. 
Exact  solutions  to  the  BGK-Boltzmann  equation  have  been  obtained  for  several 
cases  [24]  including  uniform  shear  flow  [25].  Also,  the  BGK  model  has  been 
proved  to  provide  an  accurate  representation  of  the  Boltzmann  equation  for  uni- 
form shear  flow  [26].  A  primary  new  result  of  this  thesis  is  the  development  of  a 
corresponding  kinetic  model  for  the  equations  describing  fluctuations.  The  result 
is  then  a  more  practical  coupled  set  of  equations  for  all  properties  of  interest  in 
general  nonequilibrium  states.  Here,  they  are  applied  to  the  uniform  shear  flow 
state  and  the  model  equations  are  solved  for  arbitrary  shear  rates. 

This  kinetic  description  of  fluctuations  and  transport  properties  is  at  the  de- 
tailed level  of  the  single  particle  phase  space  (position  and  velocity).  However,  it 
incorporates  as  well  a  "hydrodynamic"  description  for  large  space  and  time  scales. 
If  we  are  interested  in  space  and  time  scales  which  are  large  compared  to  the  mean 


free  path  and  mean  free  time  (distance  and  time  between  collisions),  it  is  pos- 
sible that  the  system  can  be  described  in  terms  of  few  physical  observables,  the 
local  conserved  densities  (mass,  energy  and  momentum  density).  The  averages 
of  these  local  densities  are  called  hydrodynamic  variables.  A  closed  dynamics  for 
the  hydrodynamic  variables  can  be  obtained  from  integration  of  the  Boltzmann 
equation  multiplied  by  the  corresponding  single  particle  functions  (1,  v,w2).  The 
hydrodynamic  equations  can  be  linearized  around  the  reference  state  to  deter- 
mine the  hydrodynamic  modes  (excitations)  which  describe  the  response  of  the 
state  to  the  small  perturbations.  In  Chapters  4,  5  and  6,  this  hydrodynamic 
description  is  constructed  for  states  far  from  equilibrium.  First,  the  distribution 
function  at  and  near  the  reference  state  of  uniform  shear  flow  is  calculated  for 
arbitrary  shear  rate.  This  solution  then  is  used  to  derive  a  closed  set  of  general- 
ized nonlinear  hydrodynamic  equations.  This  gives  a  complete  description  of  all 
average  properties  of  interest.  A  nonequilibrium  fluctuation-dissipation  relation 
is  used  to  express  the  dynamics  of  fluctuations  of  the  hydrodynamic  variables 
in  terms  of  the  hydrodynamic  modes.  For  mode  excitations  near  the  uniform 
shear  state,  a  long  wavelength  instability  is  found  such  that  at  any  finite  value  of 
the  shear  rate,  modes  of  sufficiently  long  wavelength  exhibit  initial  growth  rather 
than  relaxation.  Finally,  long  range  spatial  correlations  at  large  shear  rates  are 
investigated  by  calculating  the  density  fluctuations  using  these  same  hydrody- 
namic modes.  However  there  is  a  limitation  on  the  length  scale  for  studying 
spatial  correlations  because  of  the  long  wavelength  instability. 
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The  presentation  of  this  thesis  is  organized  as  follows.  In  Chapter  2,  the  exact 
formal  solution  to  the  Liouville  equation  for  uniform  shear  flow  is  derived.  To 
illustrate  an  application,  the  density  fluctuations  are  calculated  in  detail  for  a 
low  density  gas  at  small  shear  rate.  The  difficulties  of  applying  this  approach  at 
large  shear  rate  are  noted.  In  Chapter  3,  the  kinetic  theory  for  transport  and 
fluctuations  at  low  density  is  recalled,  and  three  equations  for  corresponding  ki- 
netic models  are  constructed.  The  formal  solutions  to  these  model  equations  are 
obtained.  To  illustrate  the  practicality  of  the  kinetic  model,  two  examples  are 
discussed.  In  the  remaining  chapters,  the  three  kinetic  model  equations  are  spe- 
cialized to  the  case  of  uniform  shear  flow.  In  Chapter  4,  the  stationary  solution, 
/„,  is  determined  exactly  for  all  shear  rates  and  then  a  special  "normal"  solution 
near  the  stationary  state  is  constructed.  The  normal  solution  allows  calculation 
of  the  heat  and  momentum  fluxes  in  terms  of  the  hydrodynamic  gradients  and 
identification  of  the  hydrodynamic  equations.  The  transport  coefficients  also  are 
calculated  in  detail  from  the  normal  solution.  These  fluxes  and  new  transport  co- 
efficients describe  transport  properties  far  beyond  the  usual  Navier-Stokes  order 
(e.g.  they  describe  strong  Theological  effects).  In  Chapter  5,  two  time  correlation 
functions  are  examined.  These  functions  represent  the  linear  response  to  a  small 
perturbation  of  uniform  shear  flow.  With  these  functions,  both  the  modes  of 
relaxation  and  the  stability  of  uniform  shear  flow  are  analyzed.  The  dispersion 
relations  are  solved  to  determine  these  modes  as  a  function  of  the  shear  rate.  A 
new  long  wavelength  instability  of  the  uniform  shear  flow  state  is  found.  This 
instability  is  compared  with  computer  simulation  results,  showing  an  excellent 


agreement  between  theory  and  simulation  during  the  initial  stages.  In  Chapter  6, 
equal  time  correlation  functions  are  studied.  The  long  range  spatial  correlations 
in  the  density  fluctuations  are  investigated  again,  now  at  large  shear  rates.  An 
algebraic  decay  is  expected  from  earlier  calculations  based  on  small  shear  rate, 
due  to  coupling  of  the  hydrodynamic  modes  at  equilibrium.  These  results  are 
extended  to  large  shear  rate  using  the  generalized  hydrodynamic  modes  of  Chap- 
ter 5.  The  results  are  summarized  in  Chapter  7,  and  topics  for  further  study  are 
discussed. 


CHAPTER  2 
NONEQUILIBRIUM  STATISTICAL  MECHANICS 

In  this  chapter,  a  formal  solution  to  the  Liouville  equation  with  a  noncon- 
servative  force  is  derived  for  uniform  shear  flow.  The  solution  to  this  equation 
is  used  to  calculate  fluctuations  in  this  nonequilibrium  state.  To  illustrate  its 
practical  application,  the  density  fluctuations  are  computed  at  low  density  and 
small  shear  rate. 

2.1     Formal  Solution  to  the  Liouville  Equation 

In  statistical  mechanics,  the  state  of  the  system  at  time  t  is  described  by 
the  N-body  distribution  function,  p{x\,  ■  ■  ■  ,xn,t).  Here  Xi  =  (qi,Pi)  and  q,  p 
are  the  coordinates  and  momenta  of  the  particles  in  phase  space,  respectively. 
The  time  dependence  of  the  function  is  determined  from  the  Liouville  equation. 
The  system  to  be  described  is  in  contact  with  external  reservoirs  and  the  asso- 
ciated interaction  between  the  system  and  the  reservoirs  will  be  represented  by 
nonconservative  forces.  The  Liouville  equation  with  the  external  force  is 


%  +  L  +  A 
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where 

L  =  Z  (—  ■  v<«  +  F-<(pa)  ■  vPa)  +  -  Y,  F(q°  -  «)  ■  vPt>,       (2.2) 

Q  =  l    V  m  /  *  Q?t/J 

A  =  ^VP.F«rt-  (2-3) 

The  Felt  is  an  external  force  and  F  represents  the  interaction  between  particles. 
By  solving  this  equation,  physical  observables  of  interest  and  their  fluctuations 
can  be  calculated  using  Eqs.(l.l)  and  (1.2).  Here,  we  look  for  the  solution  corre- 
sponding to  uniform  shear  flow.  This  flow  is  generated  by  Lees-Edwards  boundary 
conditions  [6].  These  conditions  are  the  usual  periodic  boundary  conditions  along 
the  x  and  z  axes,  but  with  shearing  conditions  at  y  =  ±L, 

px(y  =  -L)  =  px(y  =  L)  -  2mux(y  =  L), 

qx{y  =  -L)  =  qx(y  =  L)  -  2ux(y  =  L)t.  (2.4) 

These  are  not  the  usual  local  boundary  conditions  where  a  particle's  velocity 
accommodates  locally  the  velocity  of  the  wall.  Instead  these  are  nonlocal  relations 
between  particles  leaving  at  y  =  L  and  those  entering  at  y  =  —L.  To  motivate 
this,  it  is  useful  to  make  a  Galilean  transformation  to  the  local  Lagrangian  frame 
determined  by  the  macroscopic  flow  velocity  u*  =  CtyTj, 

Pa.i  =  Pa,i  -  mUi< 

Qa,i  =  qa,i  -  utt,  (2.5) 
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where  i,j  denotes  a  Cartesian  coordinate.  Then  it  is  seen  that  the  boundary 
conditions  (2.4)  become  simple  periodic  boundary  conditions.  This  is  particularly 
useful  for  molecular  dynamics  simulations  of  uniform  shear  flow. 

It  is  convenient  to  represent  the  Liouville  operator  in  the  mixed  set  of  vari- 
ables, {qQ,pQ}.  Then  the  new  Liouville  operator  L  becomes 

j    £  Pa,,  a         (      d         a  \  a 

Si  m  d1",i  \      dqa,  dpaiJ  dpa,i 

+  lHlP,(qa-<l0)^-,  (2.6) 

2  ait0  °P<*J 

A  =  £7F-^x<,,(iv>-  (2.7) 

a     OPa,, 

We  look  for  solutions  of  the  form 

p  =  p({pa},{qa-q/3}.*).  (2-8) 

which  expresses  translational  invariance  at  fixed  {pa}-  This  is  sufficient  to  assure 
that  the  macroscopic  state  is  indeed  that  of  uniform  shear  flow.  It  is  convenient 
to  represent  the  solution  in  terms  of  the  deviation  from  the  local  equilibrium 
distribution  function,  pi,  so  we  define  a  new  phase  function,  D, 

p(f,t)=Mf\i)eD<?'",  (2.9) 

where  T  is  a  phase  coordinate  {qa,pQ},  and  pi  is  a  local  equilibrium  distribution 
function,  pt  =  exp[-(q(t)  —  f}(t)H)].  The  q(t)  is  a  normalization  factor,  0  = 
l/kgT(t),  and  H  =  £Q  &  +  j  Ecjt/j  V{Qap),  where  V  represents  an 
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interaction  potential  between  two  particles,  and  qag  denotes  the  distance  between 
two  particles  a  and  /3.  Substitution  of  this  expression  (2.9)  into  the  Liouville 
equation  gives  the  equation  for  D(t). 

|  +  z)z)(i)  =  -A-(|  +  z)log/<  (2.10) 

=  -&-fl^((H)l-H)  +  i3LHJt.  (2.11) 


The  general  solution,  D(t),  is 

D{t)  =  e-tXD(t  =  0)  -  f*  dTe-[t-T)l 
Jo 


A  +  ^-({H),-H)-0(t)LH 


(2.12) 
The  time  dependence  in  the  brackets  on  the  right  side  occurs  entirely  through  the 
temperature  T(t).  Therefore,  at  long  times,  a  stationary  state  is  possible  if  we 
choose  the  external  force  to  hold  the  temperature  constant.  The  nonconservative 
force  can  be  fixed  as  being  proportional  to  the  velocity  relative  to  local  flow 
velocity, 

Fext(pa)  =  -Ap0,  (2.13) 

where  the  constant  A  can  be  determined  by  requiring  stationarity.  Since  the  only 
time  dependence  comes  from  temperature,  we  need  to  find  the  time  dependent 
average  energy  equation.  The  microscopic  energy  density  is  defined  by 

f  =  E  |Ur  -i»)  +  5  E  ^W*  -  q»).  (2.i4) 

The  average  energy  equation  can  be  obtained  by  the  integration  of  the  Liouville 
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equation  with  the  energy  density, 


/*"<! 


+  L  +  A)  ,5  =  0. 


(2.15) 


By  using  the  translational  invariant  property,  the  average  energy  equation  can 
be  written  as 


dt 


(e),  =  -ffly^),  +  (£F<,*«(Pa)  ■  P„)t 


(2.16) 


where  (•  ■  -)t  means  average  over  p(t)  like  in  (1.1),  and  ty  =  ZIa  P°'m°''  ^(r  ~~  la). 
Because  of  the  translational  invariance,  the  potential  part  of  momentum  flux  does 
not  contribute  to  the  energy  equation  for  uniform  shear  flow.  Now  the  stationary 
condition  is  imposed  on  Eq.(2.16  )  as  Jj(e)i  =  0.  The  proportionality  constant  A 
is  obtained  as 

a— ag&  (2.17) 


(2K) 
where  Tfi  is  a  volume  integral  of  tfj  and  K  is  a  kinetic  energy, 


r#  =  £ 


pa,pa 


(2.18) 
(2.19) 


Then  D[t)  from  (2.12)  can  be  stationary  by  taking  the  infinite  time  limit, 


Ds(f) = -  r 

Jo 


dre 


0^  -  0=*jj2LL(K  -  (K)) 


(2.20) 


where  "s"  means  stationary  and  Tu  is  total  momentum  flux,  Tu  —  Y,a  ~ 
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|  T,a^s(l0'  ~  g«)    al      •  T*"s  Ds  gives  a  formal  stationary  solution  to  the  Liou- 


2  A-a^/jWfli       laiy     a,aj 

ville  equation  for  uniform  shear  flow, 

Ps  =  PeeD-,  (2.21) 

where  pe  is  an  equilibrium  distribution  function.  By  using  this  stationary  so- 
lution, it  is  possible  to  calculate  the  average  of  the  physical  variables  and  their 
fluctuations  from  Eqs.(l.l)  and  (1.2)  in  principle.  As  an  illustration,  in  the  next 
section  we  will  discuss  the  equal  time  density  fluctuations. 

2.2     Density  Fluctuations 

In  the  last  section,  we  constructed  the  formal  solution  to  the  Liou ville  equation 
for  uniform  shear  flow.  With  this  solution,  in  this  section,  we  will  calculate  the 
equal  time  density  fluctuation.  This  fluctuation  can  be  written  by  using  Eq.(1.2) 
and  the  related  physical  variable,  density,  n  =  £Q<5(qQ  —  r).  Then  the  density 
fluctuation  can  be  represented  in  terms  of  the  correlation  function, 

(5n(T)Sn(T'))ne,  (2.22) 

where  (•  •  •)„<,  denotes  average  over  the  stationary  nonequilibrium  distribution 
function,  p„.  This  representation  is  sometimes  called  density  correlation  functions 
and  it  has  the  same  meaning  as  density  fluctuations.  These  density  correlation 
functions  can  be  rewritten  in  terms  of  the  phase  function  Ds  by  using  the  ps  = 
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peeD'  (2.21), 

(<5rc(r)5n(r')>ne  =  {5n{r)8n{v')eD')e,  (2.23) 

where  (•  •  •),  means  average  over  the  equilibrium  distribution  function,  pe,  and 
actually  it  is  the  Gibbs  distribution  function.  Now  the  density  fluctuations  are 
expressed  in  terms  of  three  phase  functions,  averaged  over  the  equilibrium  dis- 
tribution. This  representation  is  suitable  for  standard  equilibrium  molecular  dy- 
namics simulations,  because  the  equilibrium  distribution  is  well  known.  However, 
it  is  not  practical  for  theoretical  calculations  except  for  states  near  equilibrium, 
since  the  factor  eD*  is  very  complicated. 

For  states  near  equilibrium  (for  small  a),  a  perturbation  method  can  be  ap- 
plied to  the  solutions  (2.20)  and  (2.21).  Then,  the  first  term  in  (2.20)  will  be 
dominant  (the  second  term  is  proportional  to  o2),  and  the  distribution  function 
can  be  written  as 

ps  =  psl{l-0aj     dTTly(r)),  (2.24) 

where  Txy(r)  =  e~TLTxy.  For  states  near  equilibrium,  the  density  fluctuations  are 

(<5n(r)<5n(r'))„e  =  (<Sn(r)(Jn(r'))e  -  0  a  j°°  dr(Sn{r)6n(T')fxy{r))e.        (2.25) 

The  first  term  represents  an  equilibrium  density  fluctuation  and  is  usually 
written  as 

{6n{T)Sn(r'))c  =  TtS(t  -  r')  +  n2[g(r)  -  1],  (2.26) 

where  n  is  the  average  density  at  equilibrium  and  g(r)  is  the  pair  correlation 
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function.  At  low  density,  the  second  term  can  be  neglected,  so  the  density  fluc- 
tuations at  equilibrium  have  spatially  very  short  range  correlations.  The  second 
term  of  Eq.(2.25)  involves  three  phase  functions  and  is  represented  with  a  time 
integral  which  is  similar  to  the  Green-Kubo  formula.  The  Green-Kubo  formula  is 
involved  with  time  integration  of  the  flux-flux  time  dependent  correlation  func- 
tion and  gives  a  related  transport  coefficient,  while  our  second  term  has  a  time 
integral  of  two  density  fluctuations  at  one  time  and  a  flux  at  the  other  time  cor- 
relation and  it  gives  a  stationary  nonequilibrium  density  fluctuation.  We  can  see 
from  this  representation  that  the  dynamical  origin  is  very  important  for  deter- 
mining a  new  effect  of  the  stationary  nonequilibrium  fluctuation.  This  point  will 
be  more  clear  if  the  second  term  of  (2.25)  is  written  without  time  integration.  It 
is  possible  to  do  so  by  using  the  relationship 


TtJ(t)  =  ^My(t),   Mi3(t)  S  Y. <hj{t)p*&)-  (2-27) 


Then  the  density  fluctuation  for  uniform  shear  flow  is  written  as 

(«5n(r)<Mr'))„e  =  (6n(r)6n(r'))e-f3alim(6n{r)6n{r')[Mxy(t)-Mxy{C))})e.  (2.28) 

From  this  expression,  it  is  clearly  seen  that  equilibrium  "time  dependent"  fluc- 
tuations are  the  origin  of  the  effect  of  the  nonequilibrium  steady  correlation  such 
as  spatial  long  range  correlations. 
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Even  though  the  formal  representation  is  given  for  density  fluctuation  like 
(2.28),  it  is  very  hard  to  solve  analytically  because  of  the  time  dependence  which 
comes  from  the  formal  Liouville  operator.  However,  at  low  density,  the  well 
developed  kinetic  theory  is  available  to  calculate  the  density  fluctuations.  Here 
we  will  use  the  recipe  from  the  paper  by  Dufty  [27]  and  Groome,  Gubbins  and 
Dufty  [28]  to  calculate  the  density  fluctuation  without  detailed  explanations.  At 
low  density,  the  time  dependent  correlation  can  be  represented  as 

(6n(T)6n{r')Mxy{t))e  =  J  dqdp  f0pxqyTt{S{r  -  q)<5(r'  -  q')} 
-  jf'  dr  J  dqdp  foPxqyT^TnJ,[TTS{i  -  q),TT<S(r'  -  q')],  (2.29) 

where  now  Tt  depends  on  the  inhomogeneous  linearized  Boltzmann  collision  op- 
erator instead  of  a  formal  Liouville  operator, 

T,  =  exp-(v  V  +  nI)t.  (2.30) 

The  J  is  a  Boltzmann  collision  operator 

J\J,S]  =/rfvi(v  -  vO&didtf/Wfa)  -  /(»)p(«t)J.  (2.31) 

and  the  prime  denotes  the  stage  after  collisions.  The  /  is  a  linearized  Boltzmann 
operator, 


nIX 


■■  |dvi(v  -  v^bdbdpfoMlXiv)  +  X(Vl)  -  X'(v)  -  X'(Vl)},       (2.32) 
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and  /o(«i)  is  a  Maxwell-Boltzmann  distribution  function.  Also  Js  has  the  sym- 
metrized form 

J,[f,9}  =  J[f,9\  +  J[9,f}-  (2.33) 

The  calculation  of  (2.29)  will  be  simpler  if  the  arguments  in  the  Boltzmann 
collision  operator  J,  are  represented  in  Fourier  space.  Then,  after  a  few  steps 
of  calculations,  (2.29)  becomes 

(8n{r)6n{r')Mxy(t))e  = 
-(^h^f^l!^  fe'  »»v),5(-k,v)]) 

(2.34) 

where  (X,  Y)  defined  by 

(X,Y)  =  JdpfoXY,  (2.35) 

and  a(k,  v)  =  e-'("'-'kv).  After  very  long  and  tedious  calculations,  the  leading 
term  of  the  density-density  fluctuation  for  large  \r  -  r'\  represents  long  range 
correlations 

{5n(r)5n(v'))ne  ~  {5n(v)5n(v'))e  -  a  ±^X~^^^,  (2.36) 

where  T  is  a  sound  damping  constant  and  n  is  the  equilibrium  density.  The 
detailed  calculations  are  discussed  in  Appendix  A.  From  this  result,  we  can 
see  the  algebraic  decay  of  the  density  fluctuation  in  space.    The  exponent  of 
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the  spatial  decay  is  equivalent  to  the  result  obtained  from  the  phenomenological 
method  of  Lutsko  and  Dufty  [18].  Recently,  light  scattering  experiments  were 
performed  with  temperature  gradients  [4,29]  and  the  results  confirm  the  existence 
of  the  long  range  correlations. 

The  result  we  calculated  is  only  valid  at  the  low  density  limit  and  near  equi- 
librium. For  states  far  from  equilibrium,  we  should  include  the  second  and  higher 
order  terms  from  (2.20).  This  leads  to  increasingly  complex  expressions  involving 
multi-time,  multi-variable  correlation  functions.  At  present,  it  does  not  appear 
to  be  a  practical  approach.  In  the  next  chapter,  we  are  going  to  construct  the 
kinetic  "model"  equations  related  to  the  equations  of  transport  and  fluctuations 
to  find  a  way  to  circumvent  these  difficulties,  and  so  finally  to  be  able  to  calculate 
the  transport  and  fluctuations  for  states  far  from  equilibrium. 


CHAPTER  3 
KINETIC  MODELS  FOR  CORRELATIONS 


Kinetic  theory  traditionally  provides  an  equation  for  the  single  particle  re- 
duced distribution  function  (e.g.  the  Boltzmann  equation  at  low  density).  How- 
ever, kinetic  theory  can  be  extended  to  obtain  equations  for  correlation  functions 
as  well  by  applying  the  same  many-body  methods  used  to  derive  an  equation  for 
the  distribution  function.  In  the  first  section,  these  equations  for  the  distribution 
function  and  the  pair  correlation  functions  at  low  density  are  recalled.  While 
these  equations  provide  a  complete  description  of  nonequilibrium  phenomena  at 
low  density,  it  is  recognized  that  their  application  is  very  difficult  except  for  lim- 
ited cases  such  as  states  near  equilibrium.  To  allow  a  practical  study  of  states 
far  from  equilibrium,  a  well  known  kinetic  model  equation  for  the  distribution 
function  is  adopted.  In  the  second  section,  a  self-consistent  kinetic  model  for 
the  correlation  function  also  is  obtained  by  exploiting  the  exact  structural  rela- 
tionship between  the  kinetic  equation  for  the  distribution  function  and  those  for 
the  correlation  functions.  In  the  third  section,  these  model  equations  are  used 
to  describe  fluctuations  in  nonequilibrium  steady  states.  To  illustrate  the  prac- 
tical advantages  of  these  equations,  existing  results  for  correlations  at  and  near 
equilibrium  are  derived  simply  and  directly  in  the  last  section. 
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3.1     Kinetic  Equations  for  Correlations 

Kinetic  theory  is  one  of  the  starting  points  to  explain  nonequilibrium  transport 
properties  and  fluctuations.  It  describes  the  dynamics  of  the  system  in  terms  of 
reduced  distribution  functions.  If  the  interaction  potential  between  particles  is 
short  ranged,  and  the  system  is  limited  to  low  density,  then  the  equation  for  the 
single  particle  distribution  function  is  the  Boltzmann  equation.  The  same  low 
density  approximation  method  used  to  derive  the  Boltzmann  equation  can  be 
extended  to  obtain  kinetic  equations  for  the  equal  and  multi-time  phase  space 
correlation  functions  (see  Appendix  B).  Equal  and  two  time  correlation  functions 
are  the  most  important  quantities  to  study  nonequilibrium  fluctuations  because 
these  quantities  are  related  to  the  physical  observables  and  they  can  be  measured 
in  the  laboratory.  In  this  section,  the  equations  for  equal  and  two  time  correlation 
functions  are  reviewed  within  the  low  density  approximation.  [19,30,31] 

Certain  velocity  moments  of  the  distribution  function  are  the  hydrodynamic 
variables  (density,  temperature  and  velocity)  and,  similarly,  the  corresponding 
moments  of  the  phase  space  correlation  functions  are  the  hydrodynamic  fluctua- 
tions. Consequently  the  kinetic  theory  considered  below  can  be  used  to  obtain  a 
hydrodynamic  description  as  well.  This  is  discussed  in  Chapters  4-6. 

The  system  considered  here  is  a  low  density  gas  of  point  particles  interacting 
via  a  short  ranged  potential.  Let  y  =  {r,  v}  denote  the  position  and  velocity  of  a 
point  in  the  single  phase  space.  A  set  of  the  microscopic  densities  is  introduced 
as 

Hy,t)  =  J2S(y-xt(t)),  (3.1) 
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f(y,t-y',t')  =  f(y,t)f(y',t'),  (3.2) 

/(jh.Sfe,  t)  =  f(yut-  y2,  t)  -  %,  -  Vi)Htti,  t),  (3.3) 

where  Xi(t)  =  {q;(i),  Vi(t)}  denotes  the  position  and  velocity  for  particle  i  at 
time  t.  The  single  particle  phase  space  density  given  in  (3.1)  can  be  generalized 
to  a  multi-point  density  according  to 

/w(Vi,-  ■-,»..*)=     £    ■•:E*fl*(»»^4.(*))-  (3-4) 

The  average  of  this  microscopic  density  is  the  s-particle  reduced  distribution 
function, 

/Mfe, -•-»„*)■  (/Ufa,  —  .*,*»■  (3.5) 

The  angle  brackets  denote  an  ensemble  average  over  some  specified  initial  ensem- 
ble. For  simplicity,  we  chose  the  initial  time  to  be  zero  without  loss  of  generality. 
Correlations  at  two  phase  points  and  times  are  defined  by 

C(y,t;y',t')  =  f(y,f,y',t')  -  f(y,t)f(y',t') 

=  ((f(y,t)-(f(yJ)))f(y\t')).  (3.6) 

Often  a  single  time  pair  correlation,  G(yi,  3/2,  t),  is  used  for  equal  time  correlations 
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instead  of  C(y,t\y',t), 

G(yuy2,t)  =  f(yuy2,t)-  f(yi,t)f(y2,t) 

=  C(yu  t;  jft,,  *)  -  %!  -  y2)f  (yi,  t).  (3.7) 

Since  G(yi,y2,t)  is  the  correlation  function  without  the  6(yi  -  y2)  function,  it 
represents  the  "true"  correlation.  We  note  that  for  observables  of  the  form, 

A(t)«5>,   B(t)  =  '£bi,  (3.8) 

1=1  M»J 

their  averages  and  fluctuations  are  related  to  Z'1'  and  C(y,t\y',t')  by  Eqs.(l.l) 
and  (1.2)  according  to 

(A(t))t0  =  j  dya(y,t)fW{y,t),  (3.9) 

(5A(t)6B(t'))t0  =  J  dydy'a(y,  t)b(y',  t')C(y,  t;  y>,  <).  (3.10) 

The  reduced  distribution  functions  obey  the  Bogoliubov-Born-Green-Kirkwood- 
Yvon  (BBGKY)  hierarchy.  The  equations  of  the  hierarchy  can  be  closed  by  using 
a  systematic  expansion  in  powers  of  a  parameter  that  is  small  at  low  density. 
These  low  density  approximations  provide  the  equations  for  transport  and  corre- 
lation functions  at  equal  and  two  time.  The  detailed  derivation  can  be  found  in 
Appendix  B.  The  kinetic  equations  are  the  following, 


|  +  v-V 


( 
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_  +  v  .  V  -  A  )  C(y,  t;  y',  t') 


[di  +  Vl  '  Vl  "  Al  +  V2  '  V2  "  A2j  C(2/l' t;  V2' l)  =  B^'  »I/W).        (3-13) 
where  J(y\f(t))  is  the  nonlinear  Boltzmann  collision  operator 

J(y\f(t))  =  fdy'T(l,2)f(y,t)f(y',t).  (3.14) 

And  T(l,2)  is  a  binary  scattering  operator  denned  in  (B.2).  These  three  kinetic 
equations  are  a  fundamental  basis  for  the  theoretical  study  of  the  transport  and 
fluctuations  in  general  nonequilibrium  states  at  low  density.  In  principle,  these 
equations  are  solvable  because  they  are  linear  equations. 

The  Eq.(3.11)  is  the  well  known  Boltzmann  equation,  and  it  gives  a  familiar 
solution,  the  Maxwell  Boltzmann  distribution  at  equilibrium.  Since  this  does  not 
couple  to  the  other  equations,  it  is  self-solvable  if  an  initial  condition  and  proper 
boundary  condition  are  given.  However,  it  is  very  difficult  to  solve  practically 
except  near  equilibrium  because  of  the  nonlinear  collision  operator. 

The  kinetic  equation  for  the  two  time  correlation  function  in  (3.12)  applies 
for  t  >  t1  and  it  represents  a  response  at  time  t  to  a  small  perturbation  from  the 
state  which  is  described  with  f(y',t').  An  operator  A  in  (3.12)  has  a  functional 
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relationship  with  the  Boltzmann  collision  operator  such  as 

A%)  =  /^^§ff%t).  (3.15) 

From  this  relation,  we  can  see  that  the  equation  of  two  time  correlation  functions 
has  same  structure  as  that  for  the  distribution  function  linearized  around  /(f)  as 
the  reference  state.  This  also  verifies  a  generalization  of  Onsager's  assumption 
on  the  regression  of  fluctuations  in  the  nonequilibrium  state  [32,33].  In  order  to 
solve  Eq.(3.12),  an  initial  condition  should  be  provided  which  is  the  equal  time 
correlation  function,  C(yi,t;  j/2,  t).  The  equation  for  the  equal  time  correlation 
function  C{y\,t\  2/2,  t)  is  derived  from  (3.13).  This  equation  consists  of  two  of  the 
same  A  operators  and  a  new  source  term,  B(y\,  1/2,  (),  which  is  given  by 

B(j/i,y2|/)  =  7(yl,2/2l/)  +  '5(yi-!/2)/(2/i|/)-(Ai  +  A2)%I-J/2)/(yi,4),  (3.16) 

where  7(2/1,  2/2,  t)  is  defined  by 

7(2/i,2/2l/)  sr(l,2)/(yi,t)/(;/2,i),  (3.17) 

and  T(l,  2)  defines  the  binary  scattering  operator  for  a  two-particle  collision  (see 
Eq.(B.5)).  For  local  equilibrium  state,  /  =  ft,  it  is  easily  seen  from  the  Eq.(B.5) 
that  7(2/1, 2/2!/*)  vanishes  and  so  does  J(y\fi).  Then  the  remaining  source  term 
B{y\,y2\fe)  simplifies  to 

B(yi,V2\fi)  =  -(Ai  +  A2)<%i  -  riMyi,  t).  (3.18) 
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This  equation  shows  a  Fluctuation  (B)  and  Dissipation  (A  operator)  relation 
in  a  local  equilibrium  state,  in  agreement  with  the  one  derived  by  Bixon  and 
Zwanzig  [34]  and  by  Fox  and  Uhlenbeck  [35].  However  in  a  nonequilibrium  state, 
7(2/1, 2/2 1 /W)  and  <5(2/i  —  Vi)J{y\ \f{t))  are  no  longer  zero.  These  nonvanishing 
terms  are  responsible  for  different  nonequilibrium  correlation  effects.  Indeed 
7(2/1. 3ft|/(*))  is  a  true  source  term  of  unusual  nonequilibrium  fluctuations  such 
as  long  range  spatial  correlations.  This  long  range  correlation  effect  will  be  stud- 
ied further  in  section  (3.4.2).  The  meaning  of  the  true  source  term  can  be  seen 
clearly  if  the  kinetic  equation  for  G(yi,t/2,*),  defined  in  (3.7)  is  derived.  The 
correlation  G(j/i,  2/2!/)  is  a  pair  correlation  function  and  can  be  called  a  true  cor- 
relation function  since  it  is  zero  at  local  equilibrium  state.  The  kinetic  equation 
for  G(yi,y2,t)  is 


-  +  v,  •  V,  -  A,  +  v2  ■  V2  -  A2  ]  G(yu  2/2,  t)  =  7(2/1, 2/2I/M).  (3.19) 


The  inhomogeneous  term  7(2/1,2/2!/)  plays  the  role  of  a  source  for  the  correla- 
tion functions  G(i/i,2/2,i)  •  This  source  term  defined  in  (3.17)  represents  veloc- 
ity change  through  a  single  binary  collision  without  spatial  correlations.  In  a 
nonequilibrium  state  it  represents  a  violation  of  a  detailed  balance  and  generates 
long  range  nonequilibrium  fluctuations. 

Since  the  kinetic  equations  for  correlations  are  all  functionals  of  the  solution 
to  the  Boltzmann  equation,  once  the  solution  to  the  Boltzmann  equation  is  found, 
the  correlations  also  can  be  determined  in  general.  However  the  set  of  equations 
is  still  highly  complicated.  In  next  section,  we  will  use  a  BGK  kinetic  model  for 
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the  Boltzmann  collision  operator,  and  this  model  will  be  used  to  construct  the 
associated  kinetic  models  for  correlations.  These  three  model  kinetic  equations 
admit  more  practical  applications. 

3.2     Kinetic  Models  for  Correlations 

Exact  solutions  to  the  Boltzmann  equation  are  very  limited  because  of  the 
nonlinear  collision  operator  term.  Therefore  kinetic  models  have  been  proposed 
for  the  Boltzmann  collision  operator  by  replacing  it  with  a  simpler,  more  tractable, 
operator.  The  most  famous  model  of  these  is  a  single  relaxation  time  model  due 
to  Bhatnager,  Gross  and  Krook  (BGK  model)  [23].  The  essential  qualities  of 
this  model  are  its  preservation  of  the  exact  equilibrium  solution  and  all  five  con- 
servation laws.  In  this  section,  the  BGK  model  for  the  Boltzmann  equation  is 
extended  to  a  model  for  the  dynamics  of  correlations. 

When  the  kinetic  models  for  kinetic  equations  are  developed,  several  exact 
properties  of  the  kinetic  equations  (3.11)-(3.13)  should  be  preserved  in  addition  to 
their  structural  relationships.  A  first  essential  property  of  J(y\f)  is  its  stationary 
point  at  the  local  equilibrium  distribution  fi, 

J(y\ft(t))  =  0,  (3.20) 

where  fi  has  the  form 


fl(y,t)=expr£z0l{T,i)4,a(v)). 
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The  five  functions  za(r,  t)  are  arbitrary  at  this  point,  and  ^>a(v)  are  the  collisional 
invariants  of  the  nonlinear  Boltzmann  operator, 

</>c(v)  =  {l,vy}.  (3.22) 

A  second  property  of  the  collision  operator  follows  from  the  conservation  laws, 

Jdv<l>a(v)J(y\f(t))  =  0.  (3.23) 

These  two  properties  of  the  Boltzmann  collision  operator  also  are  relevant  to  the 
kinetic  equations  for  correlations  via  the  relationship  (3.15), 

AM/* A) -I,  (3-24) 

Jdvil>a{v)\{y,t\f,h)  =  0.  (3.25) 

Finally,  the  orthogonality  condition  for  the  collisional  invariants  of  the  source 
term  B(jd,yg|/(t))  is  an  important  property  to  be  preserved,  and  is  related  to 
the  conservation  laws, 

JdVlxPa(vi)B(yi,y2\f{t))  =  0.  (3.26) 

The  above  three  properties  should  be  preserved  in  constructing  kinetic  models. 

For  the  Boltzmann  collision  operator,  the  BGK  kinetic  model  preserves  the 
primary  properties  of  collision  operator,  stationarity  (3.20)  and  the  conservation 
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laws  (3.23), 

JBGK(y,\f(t))  =  -y{f(t)  -/<(*)),  (3.27) 

where  v  is  a  collision  frequency  which  depends  on  the  particle  interaction  law. 
The  conservation  laws  are  imposed  by  choosing  the  arbitrary  functions  z0(r,  ()  in 
the  local  equilibrium  distribution  to  satisfy 

fdviPa(v)U(y,t)-ft(y,t))  =  0.  (3.28) 

This  provides  five  equations  defining  the  za  as  functionals  of  /. 

A  corresponding  model  for  the  operator  A  is  obtained  directly  from  (3.15)  to 
give 

ABGK%)  iE  jdVl  Lj£-^J^{y\f[t))j  %l) 

=-"(%)-/<fei%i)jj^/iG<.*)V     (3-29) 

The  functional  dependence  of  fi  on  /  occurs  entirely  through  the  za.  Therefore, 

-Jt{y,  t)  =  fi(y,  t)M^jJ—rMT'  *)■  (3-30) 


*/(ih,0  ;<5/(yi,«) 

A  summation  convention  over  repeated  Greek  labels  is  assumed.  The  functional 
derivative  in  (3.30)  is  evaluated  by  using  the  conservation  law  property  (3.28) 
and  it  gives  the  BGK  model  for  A  operator 

ABGKh(y)  =  -v(l-P)h(y),  (3.31) 
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where  P(t)  is  a  projection  operator  defined  over  functions  X(y)  by 

'(»)  =  My,  t)i>a(y)g  $(',  *)  /  *i  to(vi)2T(r,  vx).  (3.32) 


PXf 


The  detailed  derivation  for  ABGK  is  in  Appendix  C. 

Next,  for  the  equal  time  correlation  function,  a  model  for  the  source  term 
B(y\,V2\f)  is  required  in  addition  to  ABGK.  The  primary  constraint  for  con- 
structing _gBGK  is  the  orthogonality  condition  in  (3.26).  Therefore  we  first  make 
the  orthogonality  of  B  to  the  collisional  invariants  explicit  by  inserting  two  or- 
thogonal projections  in  (3.16), 

B(vi,  »l/)  =  (i  -  Pi)(i  -  ft)W*»,»l/J  +  i(tn  -  ¥i)Jinn\f) 

-(A1+A2)S(yl-y2)f(yi)\-  (3.33) 

This  result  is  still  exact.  To  define  the  corresponding  BGK  model  for  B,  we 
replace  J  and  A  by  JBGK  and  ABGK  and  neglect  (1  -  Pi)(l  -  P2)7(j/i,  Vi\f)  (it  is 
a  choice  we  make  for  simplicity),  leading  to  the  result 

BB0K(IA,I&I/)  =  "(1  -  ft)(l  "  PiWvi  ~  M)(/te)  +  fl(Vl)),  (3-34) 

where  use  has  been  made  of  the  property  (l-P1)(l-P2)(Pi  +  P2)  =  0.  While  the 
orthogonality  condition  is  not  sufficient  to  determine  BBGK  uniquely,  the  choice 
(3.34)  is  consistent  with  the  known  form  near  equilibrium  (3.18).  With  this  BBGK 
model,  the  kinetic  model  for  the  source  term  7BGK  in  the  equation  for  G  (3.19) 
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can  be  defined.  The  definition  for  B  (3.16)  is  used 

7BGK(yi,2/2|/)  =  BBGK(yi,y2\f)  -%,  -ift)JBGK(»,|/) 

+(AfGK  +  A2BGK)<5(y1-!/2)/(yi).  (3.35) 

Then  the  model  for  7BGK  is 

7BGK(yi,y2|/)  =  vHM»  -  V2){f(yi)  -  Mvi)},  (3.36) 

where  use  has  been  made  of  the  property  Pi(l  -  fy)ft{Vi)S(Vi  -  Vi)  —  0  ■ 

Now  the  extended  BGK  model  for  kinetic  equations  has  been  completed. 
Eqs.(3.11)-(3.13)  are  now  replaced  by 

{e + v ' v)  /(y' ()  =  ~v  (/(y' t]  ~  /,(y' t)] '  (3'37) 

(^  +  C(t))c(y,t;y',t')  =  0,  (3.38) 

(^  +  C1(t)  +  C2(t)\c(y1,t;y2,t)  =  u(l-Pl)(l-P2)5(yi-y2)(f(yl,t)+fl(yut)l 

(3.39) 
where  the  operator  C  is  defined  by 

£(i)  =  v-V  +  j/(l-P).  (3.40) 
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The  equation  for  G(y1,y2,t)  also  is  written 


JUAM  +  AW     3(J/l,t;y2,t) 


Equations  (3.37)-(3.39),  or  (3.41)  provide  the  transport  properties  and  fluc- 
tuations in  nonequilibrium  gas.  Since  these  equations  are  derived  without  any 
restriction  to  near  equilibrium  states,  these  equations  can  be  used  for  practical 
applications  of  transport  and  fluctuations  far  from  equilibrium. 

The  equations  can  be  solved  by  first  solving  the  BGK-Boltzmann  equation 
for  given  initial  and  boundary  conditions.  Once  the  solution  to  the  transport 
equation  is  found,  it  is  used  to  determine  the  operator  C  and  the  source  term  B 
or  7.  To  solve  the  two  time  correlation  function  equation  for  t  >  t',  for  the  initial 
value  is  needed  at  t  =  t'.  This  initial  value  is  calculated  from  the  solution  to 
(3.39)  or  (3.41).  In  the  next  section,  the  formal  solutions  to  the  equations  for  the 
correlations  are  constructed  with  the  assumption  that  the  solution  to  the  BGK- 
Boltzmann  equation  is  known.  The  solution  to  the  BGK-Boltzmann  equation 
will  be  calculated  in  Chapter  4  for  uniform  shear  flow  as  a  special  nonequilibrium 
state. 

3.3     Steady  State  Correlations 

The  kinetic  models  developed  above  allow  the  study  of  correlations,  once  the 
solution  to  the  Boltzmann  equation  is  obtained.  BGK-Boltzmann  equation  is 
still  a  highly  nonlinear  integro-differential  equation  due  to  the  dependence  of 
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ft  on  /.  However,  exact  solutions  of  a  Boltzmann  equation  for  nonequilibrium 
steady  states  corresponding  to  heat  and  momentum  transport  have  been  obtained 
recently  for  planar  geometries  [36].  Since  the  system  considered  in  this  thesis 
is  uniform  shear  flow,  the  exact  solution  to  BGK-Boltzmann  equation  will  be 
studied  specifically  in  the  next  chapter.  In  this  section,  to  illustrate  the  formal 
solutions  to  the  correlation  function  equations  (3.38,  3.39),  we  assume  that  a 
steady  state  distribution,  /„  is  known.  Then  the  formal  solution  to  Eq.(3.38)  is 

Cs(y,  t;  y',  0)  =  J  dys  K(y,  Jft, t)C.(jft,  0;  y',  0),  (3.42) 

K(y,y',t)  =  exp  {-Ct}  S(y  -  y').  (3.43) 

Stationarity  has  been  used  to  choose  t'  =  0.  The  initial  condition  in  (3.42)  is 
determined  from  the  stationary  solution  to  the  equal  time  correlation  function 
equation,  (3.39), 

Cs(y,y')  =  Cs(y,0\y\0) 

=  J™  dr  J  dy1dy2K{y,y1,T)K(y',y2,T)BBGK(yuy2\fs).    (3.44) 

The  two  formal  solutions  to  the  correlation  equations  are  expressed  in  terms  of 
the  Green's  function  K(y,y',t),  defined  in  (3.43).  The  equation  for  this  Green's 
function  is 

(it  +  V  '  V  +  ")  K{V' V'' 4)  =  "AtvM.M^Of.  y'. f),  (3.45) 
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with  the  initial  condition  K(y,y',0)  =  S(y  -  y').   The  functions  Ka(r,y',t)  are 
defined  by 

Ka(r,y',t)  ««5ftr|/.)y*riMv,)Jir0r,Y1,¥',*).  (3.46) 

The  solution  to  (3.45)  can  be  reduced  to  the  corresponding  ideal  gas  problem  as 
follows.  Integration  of  (3.45)  gives 

K(y,y',t)  =  K0(y,y',t) 

-v  J  dr  J  dyi  K0(y, yx,t  -  ^f^y^ip^v^K^T^,  y', r), 

(3.47) 


+  1 

Jo 


where  we  have  introduced 


K0(y,  y',t)  =  exp  [-  (v-  V  +  v)  t]  5{y  -  y').  (3.48) 

The  remaining  set  of  five  functions  Ka(r,  y1,  t)  is  determined  by  substituting  (3.47) 
into  the  right  side  of  (3.46)  to  give 


Kc 


(r,  y\  t)  =  K0a{i,  y',  t)  +  J  dr  J  drJag(T,  tut-  T)Kff(iuy',  r),       (3.49) 
(r,j/,t)  =  g$(t\f.)  J 'dv^(v)*b(tf,*',«),  (3.50) 


with 
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and 

Iaff(v,  n,t)sv £j(r|/,)  J  dvdv,  iP-,(v)K0(y,  yu  t)fs(Vl)M^)-  (3-51) 

The  quantity  K0  is  determined  entirely  from  the  dynamics  of  the  free  streaming 
gas  .  Then  KQa  and  Iag  are  obtained  from  integration,  and  the  Green's  function  K 
is  also  calculated  in  (3.47).  Finally,  the  equal  time  correlation  function  (3.42)  and 
dynamic  correlation  function  (3.44)  can  be  calculated.  Examples  of  this  analysis 
are  shown  in  the  following  sections  for  equilibrium  states  and  stationary  states 
near  equilibrium.  From  the  confirming  results  of  these  examples,  the  validity  of 
kinetic  models  for  BBGK  and  7BGK  is  supported. 

3.4     Examples 

Two  simple  states  are  considered  to  illustrate  the  analysis  of  kinetic  models 
in  previous  section.  First  the  equilibrium  state  is  considered,  and  from  the  same 
systematic  analysis  we  calculate  the  extended  hydrodynamic  modes.  The  second 
example  concerns  a  stationary  state  near  equilibrium.  The  perturbation  method 
about  hydrodynamic  field  gradients  is  used  to  get  a  stationary  state,  and  the  equal 
time  correlation  function  is  calculated.  In  particular,  the  results  of  Chapter  2  for 
long  range  spatial  correlations  are  recovered. 
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3.4.1     Equilibrium  States 

At  equilibrium  the  distribution  function  is  the  local  equilibrium  state  distri- 
bution. It  is  actually  the  Maxwell-Boltzmann  distribution  function, 

(771        \  ^^ 
^-^J      exp(-mv2/2kBT),  (3.52) 

where  n  is  the  density  m  is  the  mass  of  the  particle,  ks  is  the  Boltzmann  constant, 
and  T  is  the  temperature.  The  source  term  BBGK  given  by  (3.34)  simplifies  to 

BBGK(yi,2/2|/e)  =  -(A  +  C2)8(yi  -  a,)/.**,),  (3.53) 

in  agreement  with  (3.18).    Direct  calculation  from  (3.43)  and  (3.44)  gives  the 
expected  results  for  the  low  density  equilibrium  correlations, 

Ce(y;y')  =  6(y-y')My).  (3.54) 

Use  of  this  as  an  initial  value  to  (3.42)  gives  the  two  time  correlation  functions 
as 

Ct(y,t;tf,0)  =  KM,t)f.(tf).  (3.55) 

The  Green's  function  is  determined  from  (3.45)-(3.49).  Fourier  -  Laplace  trans- 
formation is  used,  giving 

Ce(k,k',z;v,v')  =  I  drdr'ei(kr+k  r)  /     dte-ztCe(y,t;y' ,0) 

=  (27r)3,5(k  +  k'JC^k,  z|v,  v').  (3.56) 


38 
The  correlation  function  is  now  written  as 

Ce(k,z|v,v')  =  tf(k,Z|v,v')/e(v'),  (3.57) 

and  the  Green's  function  A"(k,z|v,v')  is 

if(k,  z|v,  v')  =  i?(k,  z\v)6(v  -  v') 

+  -R(k,  z|v)fl(k,  z\v')fe(v)<l>a(v)Mv')B^(k,  z),     (3.58) 
n 

where 

R(Mv)  =  (z  +  i/-ik-v)-\  (3.59) 

Baff{k,  z)  =  SaP  -  Ia/3(k,  z).  (3.60) 

The  functions  <f>a  are  linear  combinations  of  the  summational  invariants,  il>a  and 
are  defined  in  Appendix  D.  The  first  term  of  the  Green's  function  (3.58)  describes 
the  free  streaming  motion  with  a  pure  exponential  damping  by  v.  The  second 
term  shows  the  dynamics  of  correlations  determined  by  the  singularity  of  the 
matrix  B~^(k,  z).  The  detailed  form  of  matrix  B  is  shown  in  Appendix  D.  The 
determinant  of  the  matrix  B  gives  five  simple  poles  at  za(k)  which  originated 
at  z(k  =  0)  =  0.  A  pole  at  the  origin  represents  a  conserved  quantity,  so  the 
five  simple  poles  at  finite  k  characterize  five  hydrodynamic  modes  which  are  two 
sound  modes,  a  heat  mode  and  two  shear  modes.  Since  there  is  no  restriction 
on  k,  it  is  possible  to  calculate  these  poles  not  only  at  small  k  (Navier-Stokes 
order)  but  also  at  larger  k.  This  defines  the  extended  hydrodynamic  modes.  The 
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poles  are  calculated  numerically  for  finite  k,  and  it  is  found  that  for  small  k  the 
dynamics  of  correlations  is  equivalent  to  that  from  the  linearized  Navier-Stokes 
equations.  For  large  k,  these  extended  modes  exist  up  to  some  maximum  values  of 
k  corresponding  to  wavelengths  as  compared  to  the  mean  free  path.  The  real  and 
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Figure  3.1:  Real  and  imaginary  parts  of  the  extended  hydrodynamic  modes  from 
the  BGK  model  (in  units  of  v).  ±  indicates  sound  modes,  H  the  heat  mode  and 
S  the  shear  modes. 
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imaginary  parts  for  the  two  propagating  sound  modes  (±),  ,  the  heat  diffusion 
mode  (H),  and  the  shear  diffusion  modes  (s)  are  shown  in  Fig.3.1.  In  addition 
to  the  hydrodynamic  modes,  there  are  poles  and  branch  points  which  approach 
—v  as  k  — >  0.  These  represent  the  kinetic  modes  with  lifetimes  of  the  order  of 
i/'1  for  all  k. 

This  equilibrium  dynamic  correlations  are  well-studied  field  from  the  lin- 
earized Boltzmann  equation  using  more  complex  models  [37],  The  analysis  with 
our  kinetic  models  for  correlations  shows  that  all  qualitative  features  of  fluctua- 
tions in  equilibrium  states  are  preserved. 

3.4.2     Stationary  States  Near  Equilibrium 

As  a  second  example,  a  stationary  state  near  local  equilibrium  is  considered. 
First,  in  order  to  solve  the  BGK-Boltzmann  equation,  the  Chapman-Enskog  ap- 
proximation, which  is  a  power  series  expansion  in  the  gradients  of  the  hydrody- 
namic variables,  is  used.  The  stationary  solution,  obtained  by  assuming  small 
temperature  and  flow  field  gradients,  is  given  by 


/.(y)  -» Mv) 


-i^fsM+f^M) 


(3.61) 


Here  X(r)  is  the  temperature,  u(r)  is  the  flow  field,  and 


S(v)  ee  V(\mV*  -  -2kBT),  (3.62) 
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Dtf(v)Hro(W$-|v%),  (3.63) 

where  V(r)  =  v  —  u(r)  is  the  peculiar  velocity. 

The  equal  time  correlation  function  is  calculated  in  terms  of  an  integral  repre- 
sentation for  the  solution  to  (3.41)  rather  than  (3.39)  in  order  to  see  the  difference 
from  equilibrium  correlations  more  directly.  It  is 

Gs{yuV2)=  I"  dte-^+c^(yx,y2\}s).  (3.64) 

JO 

This  shows  a  signature  of  long  ranged  spatial  correlations  which  arise  from  the 
hydrodynamic  spectrum  of  C.  It  can  be  calculated  further  by  using  the  hydrody- 
namic  eigenfunctions  of  C  to  the  lowest  order  in  the  gradients.  The  eigenfunctions 
are  proportional  to  the  summational  invariants,  il>a,  so 


G, 


rOO 

/     die-^+^r,  -  T2)Mal,Mvi)Mv2)Mvi)MV2) 
Jo 

=  (£,  +  C2)-l6(r].  -  r2)Ma0Mvi)Mv2)fi(yi)fi{y2),       (3.65) 


where 


McfSivt  -r2)  =  gZlggl  J  dv1dv2i/-(,(v1)^(v2)7(y1,j/2|/).  (3.66) 

In  Fourier  representation,  the  eigenvalues  of  £  lead  to  terms  oc  (A„(fc)  +  Xp(k))~l 
which  lead  to  slow  algebraic  decay  in  configuration  space.  This  is  a  very  interest- 
ing phenomena  because  there  are  no  long  ranged  correlations  in  the  equilibrium 
state.  These  long  ranged  correlations  near  equilibrium  have  been  verified 
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perturbatively  [2]  and  recently  accurate  experimental  results  [4,29,38]  are  sup- 
porting the  existence  of  the  long  ranged  correlations.  The  coefficients  Mag  can 
be  calculated  easily  and  explicitly  .  First,  7(2/1,2/2!/)  is  found  from  (3.36), 


c(yuV2\!s)  =  -PiP2^yi-y2)Si(vi) 


Br^S^+fe^WJ 


Next,  the  result  for  Mag  is  calculated  to  be 
Mal)-- 


1    /r-i^c        ,   9u'n        ) 


where 


See,,  s  g~^  J dVf<(v)Mv)Mv)S>(v), 
DaP.i,  =  9^911  fdvfl(v)1>.(v)My)Dii{v)). 


(3.67) 

(3.68) 

(3.69) 
(3.70) 


These  results  are  exactly  the  same  as  those  obtained  from  the  true  kinetic  theory 
(3.19)  [39], 

With  this  general  result,  the  density  correlation  function  can  be  calculated 
for  uniform  shear  flow  in  order  to  compare  the  result  with  the  one  we  got  from 
the  Liouville  equation  in  Chapter  2.  The  density  correlation  function  is  defined 
by 

(6n(ii)5n{T2)) 
=  J dvtdViC.fa.Vt)  (3-71) 

=  n5(Tl-n)  +  Jo°°dtJdvldv2e-^+c>)t1(y1,y2\fs).  (3.72) 
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The  first  term  of  Eq.(3.72)  represents  equilibrium  density  fluctuations  and  is 
same  as  (2.26)  with  very  short  range  correlation.  The  second  term  of  Eq.(3.72) 
represents  the  second  term  of  Eq.(2.25)  at  low  density.  By  using  the  coefficients 
Map  as  defined  in  (3.68),  the  second  term  is  rewritten  as 

J™  dt  J  dVldv2  e-^+c'^(yuy2\fs) 

=  a(n  -  r2)Ma0  J™  dt  J  t*v,  e-**S»(vt)A(yi) 

xJdv2e~c>tTP(V2)ft(y2).  (3.73) 

The  coefficients  Map  in  (3.68)  have  only  one  velocity  gradient  term  in  uniform 
shear  flow, 

Ma0  =  -T^B./ifl.  (3.74) 

A  Fourier  transformation  is  performed  for  the  density  correlations.  The  oper- 
ator C  becomes  the  C  =  -ik  ■  v  +  v(l  -  P).  The  hydrodynamic  eigenvalues 
and  eigenvectors  of  the  operator  C  are  equal  to  the  hydrodynamic  spectrum  of 
the  inhomogeneous  linearized  Boltzmann  operator  which  can  be  found  in  (A.5)- 
(A.9).  Also  for  simplicity,  the  collisional  invariants  i/j„  are  chosen  to  be  the  same 
orthonormal  set  {<j>a}  as  was  used  in  equilibrium  analysis  in  the  first  example, 
Eq.  (D.2).  By  using  the  lowest  order  eigenvectors  of  C  and  the  nonzero  elements 
of  Mcff,  the  density  correlation  function  is  found  to  be 


(<5n(r1)(5n(r2)) 
=  irffc  -  r2)  -  ^-32najdke-'k^-^^ 


3 

+ 

3         TA:2        ' 

2or/t2 

io  r2fc4  +  fcv 
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(3.75) 
<■/  x       3  an  (xi  —  x^Htyi  —  V2) 

" n<5(ri  - F2)  -  80^f        (ri      r2y         +  hlgher  °rder-  (3-76) 

This  result  is  exactly  same  as  the  one  in  (2.36)  which  was  obtained  in  Chapter 
2.  This  proves  that  the  BGK-kinetic  models  preserve  a  correct  description  of 
fluctuation  phenomena  in  stationary  state  and  near  equilibrium,  and  they  also 
provide  a  simpler  and  more  practical  approach  to  the  problem. 


CHAPTER  4 
TRANSPORT  EQUATIONS  FOR  UNIFORM  SHEAR  FLOW 

Kinetic  model  equations  for  transport  and  fluctuations  were  developed  in  the 
last  chapter.  We  now  specialize  those  equations  to  uniform  shear  flow.  In  this 
chapter,  the  BGK-Boltzmann  equation  for  the  distribution  function,  /,  is  exam- 
ined to  study  transport  properties  of  the  system.  The  stationary  reference  state 
solution  will  be  determined  and  discussed.  Later,  a  "normal"  solution  for  states 
near  this  stationary  state  is  obtained  using  a  modified  Chapman-Enskog  method. 
This  normal  solution  then  is  used  to  calculate  the  average  heat  and  momentum 
fluxes,  which  obtains  a  hydrodynamic  description.  These  hydrodynamic  equa- 
tions are  accurate  to  second  order  in  spatial  gradients  relative  to  uniform  shear 
flow  state,  but  there  is  no  restriction  on  value  of  the  shear  rate.  The  shear  rate 
dependent  generalized  transport  coefficients  are  given  explicitly. 

4.1     Stationary  State  Solution 

In  the  last  chapter,  a  more  practical  and  theoretical  description  of  nonequi- 
librium  phenomena  has  been  given  by  the  BGK  model  equations  (3.37)-(3.41). 
However  the  first  of  these  equations,  the  BGK-Boltzmann  equation,  is  highly  non- 
linear and  difficult  to  solve  in  general.  One  of  the  few  cases  for  which  an  exact 
solution  is  known  is  that  for  uniform  shear  flow.  Consequently,  in  the  remainder 
of  this  thesis,  attention  will  be  restricted  to  uniform  shear  flow  and  states  near 
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uniform  shear  flow.  The  uniform  shear  state  is  a  planar  flow  whose  x-component 
of  the  flow  velocity  has  a  gradient  along  the  y-axis,  U;  =  a^rj,  a;j  =  a5ix5]y, 
where  a  is  a  constant  shear  rate.  In  addition,  the  density  n  and  temperature 
T  are  spatially  uniform.  This  state  is  generated  by  a  periodic  boundary  condi- 
tion in  the  local  Lagrangian  frame  [6],  so  there  is  no  boundary  layer,  as  there 
is  for  planar  Couette  flow  [40].  Due  to  the  simplicity  of  uniform  shear  flow 
state  at  the  macroscopic  level,  it  has  been  studied  extensively  as  a  prototype  of 
nonequilibrium  states  far  from  equilibrium  in  theory  [25,41,42]  and  in  computer 
simulations  [7-9,26].  The  boundary  driven  shear  leads  to  viscous  heating.  For 
practical  purposes,  it  is  useful  to  introduce  a  uniform  "thermostat"  to  control  this 
heating  and  establish  a  steady  state.  This  is  done  in  most  simulation  studies. 

Consider  first  the  BGK-Boltzmann  equation  for  a  general  nonequilibrium 
state, 

i  +  v  •  Vr  +  Vvm~l  •  O  /  =  -K/  -  /*),  (4.1) 

where  Fext  which  serves  as  a  thermostat,  is  an  external  force  applied  on  each 
particle.  There  are  several  thermostats  that  have  been  used  in  both  theory  and 
computer  simulations.  We  choose  here  a  force  proportional  to  the  velocity  relative 
to  the  actual  flow  field  and  it  is  in  the  form  of 

Fei((v,  t)  =  -mA(v  -  u(r,  t)).  (4.2) 

The  proportional  "constant"  A  is  determined  by  requiring  stationarity  of  the  sys- 
tem, and  may  depend  on  the  density,  temperature  and  shear  rate.  The  parameter 
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v  in  the  (4.1)  is  a  collision  frequency  which  depends  on  an  interaction  law.  This 
frequency  is  a  function  of  the  density  and  temperature.  At  low  density,  it  can  be 
written  as 

v  ~  nTb,  (4.3) 

when  the  potential  has  a  form  of  V(r)  ~  r~',  with  6=1/2-  2/1.  In  the  case  of 
Maxwell  molecules  (I  =  4),  6  is  zero  so  u  becomes  independent  of  temperature. 
For  the  hard  sphere  case,  /  tends  to  infinity  resulting  in  a  value  for  b  of  1/2.  Also, 
Si  in  (4.1)  is  the  local  equilibrium  distribution  function  of  (3.21),  written  more 
explicitly  as 

^^^^^^^j^^l-^ib^-^012)'  (4-4) 

The  hydrodynamic  fields  are  defined  by 

n(r,t)  =  |dv/(r,v,t),  (4.5) 

n(r,t)u(r,i)  =  |rfvv/(r,v,(),  (4.6) 

t(r,  t)kBT(T,  t)  =  Jdv  l-m(y  -  u(r,  t))2f(r,  v,  t),  (4.7) 


3 

2m: 


which  identifies  n  as  the  density,  T  as  the  temperature  and  u  as  the  flow  velocity. 
The  scales  of  interest  are  for  the  distances  larger  than  the  mean  free  path  and 
for  times  longer  than  the  mean  free  time.  In  this  regime,  the  system  is  charac- 
terized by  conserved  hydrodynamic  variables  (such  as  density,  momentum  and 
energy).  The  local  conservation  laws  are  obtained  from  the  Boltzmann  equation 
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by  integration  over  the  velocities  with  summational  invariants  (l,v,  v2)  to  get 

3tl 

—  +  V-nu  =  0,  (4.8) 

at 

^  +  u.Vu,  +  i^=0,  (4.10) 

dt  p  Or j 

where  n,  T  and  u  are  the  density,  temperature  and  flow  velocity.  Also,  ty  is  the 

pressure  tensor, 

Uj  =  I  dvm(vi  -  Ui)(vj  -  Uj)J,  (4.11) 

and  q  is  the  heat  flux, 

q=   f  dv)-m{v-u)2(v-u)S.  (4.12) 

The  inhomogeneous  term  w  in  the  temperature  equation  is  due  to  the  external 
nonconservative  force  representing  the  thermostat,  and  is  written  as 

w  =  -3\nkBT.  (4.13) 

The  external  force  does  not  appear  in  the  velocity  equation  since  the  average  of 
(4.2)  vanishes. 

Now  uniform  shear  flow  is  considered.  It  is  useful  to  express  the  kinetic 
equation  (4.1)  in  terms  of  the  velocity  in  the  local  rest  frame,  defined  by  v[  = 
V{  —  dijTj.  In  this  frame,  the  macroscopic  uniform  shear  state  becomes  spatially 
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homogeneous.  Consequently,  we  look  for  a  stationary  solution  to  Eq.(4.1)  of  the 
form 

/(r,v,i)  =  /5(v').  (4.14) 

The  equation  for  /,  becomes, 

where  s  denotes  stationary  and  f3i  is  a  stationary  local  equilibrium  distribution 
function  (Maxwell-Boltzmann  distribution)  with  the  hydrodynamic  fields  for  uni- 
form shear  flow.  The  solution  to  Eq.(4.15)  is 

/,(v')  =  "/     dre-UTeTLft(V).  (4.16) 

Jo 

where  L  is  the  velocity  operator, 

L  =  a^  +  Xk[-  (4'17) 

For  an  arbitrary  function  of  v,  X(v),  it  has  the  property, 

eiLX(v)  =  e3A(X(eaAy(-t)uj),  (4.18) 

where  Ay  (4)  =  <Sy  -  ay  4.  Therefore  the  stationary  solution  can  be  written 

/s(v')  =  v  f°°dt  e-^-^Me^A.A-tWi)-  (4-19) 

Jo 
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The  parameter  A  is  determined  from  the  temperature  equation  (4.9)  by  requiring 
5  =  0  and  it  is  found  that 

at 

A  =  £*%-.  (4.20) 

The  components  of  the  momentum  flux  can  be  calculated  from  Eq.(4.11)  by  direct 
integration  to  get 

t,,IX  =  ps  +  i^^Ps,  (4.21a) 

2A 
4-  =  Ps-2AT^'  (4-21C) 

is,z,  =  «.«,  (4-21d) 


where  ps  =  n„kBTs,  tSiij  =  i3ji  and  all  other  elements  of  t,#  are  zero.  Use  of  tStXy 
in  (4.20)  gives  the  equation  to  fix  X(a), 

3\(2\  +  is)2  =  ua2.  (4.22) 

This  equation  has  one  real  solution  and  two  complex  conjugate  solutions.   The 
physically  relevant  real  value  is 


A  =  ?sinh[cosh-1(l  +  9^)/6]2.  (4.23) 


From  these  expressions  for  the  momentum  fluxes,  the  shear  viscosity  17(a)  and 
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viscometric  functions  $1(0),  ^2(1)  are  determined  to  be 


^W-^— 3(£h^  (4-25) 

»,(«)  =  ^^  =  0,  (4.26) 


where  all  of  the  shear  rate  dependence  occurs  through  A  =  A(a)  in  (4.23).  The 
same  integration  can  be  done  for  heat  flux  in  (4.12),  and  it  is  found  to  vanish, 
as  expected.  Therefore  (4.24)-(4.26)  describe  all  the  transport  properties  for  the 
reference  state  of  the  uniform  shear  flow.  The  Figure  4.1  shows  rj'(a')  =  vr\lys 
and  ^l(a')  =  i/2}ii/p3  as  a  function  of  a*  =  a/ v.  The  results  agree  with  previous 
studies  [43,44].  They  exhibit  strong  rheological  effects  at  large  shear  rates,  even 
though  we  consider  a  simple  gas.  The  transport  coefficients  behave  for  small 
shear  rate  as 

/        4a2      28a4  \  ps  , 

\       3  v2       9 1/4  )  v  v       ' 

«-f-**-'&*'")*  <«' 
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Figure  4.1:  Shear  viscosity  ?7(a)(solid  line)  and  viscometric  function  *(a)  (dashed 
curve)  are  shown  as  function  of  a.  r)  is  in  unit  of  pjv,  *  is  in  unit  of  p,/v2  and 
a  is  in  unit  of  v. 


and  for  large  a  as 


<a\-^  ff3\2/3 


-2/3  1    /3\V3    /as-4/3 


*1  = 


3-3(V)      (-)       +(3)      (-)       +  ■ 


(4.29) 


.(4.30) 
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Now  all  hydrodynamic  variables  and  fluxes  are  known  for  the  stationary  uni- 
form shear  state.  Our  next  interests  are  to  study  fluctuations  in  this  state,  and 
to  see  the  behavior  of  states  which  deviate  from  the  stationary  state  by  small 
spatial  gradients.  In  the  next  section,  a  variant  of  the  Chapman-Enskog  expan- 
sion method  is  introduced  to  the  BGK-Boltzmann  equation  to  find  a  "normal" 
solution  for  these  latter  states. 

4.2     Normal  Solution 

Our  objective  here  is  to  find  a  solution  which  describes  small  perturbations 
relative  to  stationary  uniform  shear  flow  due  to  spatial  gradients  in  the  hydrody- 
namics fields.  Since  the  stationary  state  is  homogeneous  in  the  local  rest  frame, 
it  is  convenient  to  continue  using  the  relative  velocity  v'  =  v  —  us.  The  BGK- 
Boltzmann  equation  is  then  written  as 

(|  +  W  +  <w)  J:  -  JM  ■)  + » sm-^J  /(r> v'.  t) 

=  -u(f(Ty,t)-fe(r,v',t))    ,  (4.31) 

where  <5u(r,  t)  is  the  deviation  of  the  flow  velocity  from  that  of  the  reference 
state,  <5u(r,  t)  =  u(r,£)  -  u,(r).  The  operator  L(v',a)  is  same  as  in  (4.17).  We 
look  for  solutions  for  which  all  of  the  space  and  time  dependence  occurs  through 
the  hydrodynamic  fields,  ya(r,t)  =  {n(r,t),T(r,  £),<5u(r,  t)}.  Such  a  solution 
is  called  a  "normal"  solution.    Near  the  stationary  state,  the  spatial  gradients 
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of  the  hydrodynamic  variables  are  small.  Since  the  distribution  function  is  a 
function  of  the  hydrodynamic  variables  (i.e.  normal  solution),  we  can  expand  the 
distribution  function  with  the  small  gradients  as  an  expansion  parameter, 

/(yo(r,i),v')  =  /<0»(yo(r,i),v')+£/(1»(ya(r,i),v')  +  ---,  (4.32) 

where  e  is  a  formal  parameter  representing  each  spatial  gradient  relative  to  the 
reference  state.  This  is  analogous  to  the  Chapman-Enskog  method  for  states 
near  local  equilibrium,  but  here  the  reference  state  can  be  far  from  equilibrium. 
In  order  to  determine  each  /'r)  in  Eq.  (4.31),  we  need  to  represent  the  time 
derivatives  in  terms  of  spatial  gradients  as  well  using  the  hydrodynamic  equations. 
Hence  the  time  derivative  also  can  be  expanded  in  powers  of  the  same  e, 

9     9i0)     a(1)  „,„N 

By  putting  the  expansion  (4.32)  and  (4.33)  into  Eq.(4.31),  equations  for  each  /'r' 
can  be  determined.  The  zeroth  order  equation  is 

3(0)  \ 

—  -  L(v\  a)  +  A<5u  ■  V'v  1  /<">  =  -K/(0)  -  A),  (4.34) 

and  the  first  order  equation  is 

5? + «+ **>&)  /(0)  =  -  (?  -  m  «•) + ^i, + A  /»•  iu*> 
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First,  the  zeroth  order  equation  is  considered.  This  equation  has  the  same 
form  as  for  /,,  (4.15),  except  that  the  hydrodynamic  fields  can  have  a  local  space 
and  time  dependent  parts.  Therefore  the  solution  to  Eq.(4.34)  can  be  identified 
as 

exp  {-2kBT{Trf2X(''t)Tl^-TM  ~  <^M))i2)  •  (4-36) 

For  consistency,  the  hydrodynamic  variables  in  /'"'  must  satisfy  the  equations 

d<°> 

— n  =  0,  (4.37) 

%f  =  0  (4.38) 

^p  +  ai;<SUj=0.  (4.39) 

These  conditions  can  be  verified  easily  since  they  are  nothing  but  the  hydrody- 
namic equations  (4.8)-(4.10)  at  lowest  order  e, 

g(0) 

-g^n  =  0,  (4.40) 

-nAs— +  a,Ji|°)  +  3Ap  =  0,  (4.41) 

-^-i  +  aijSUj  =  0,  (4.42) 
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where  p  denotes  pressure  and  t ^    is  lowest  order  momentum  flux  defined  by 

4°>  =  Jdvm(v,  -  u,)(v3  -  U])fm.  (4.43) 

The  extra  terms  a^tfj  +  3Ap  in  Eq.(4.41)  vanish  if  we  calculate  the  values  of  tjj' 
from  (4.43).  The  (s?  are  calculated  with  a  similar  integration  method  used  in 
(4.21a)  and  they  are  the  same  as  tStij  except  that  now  the  density  and  temperature 
have  local  space  and  time  dependent  parts, 

t^=p(r,t)  +  j^p(r,t),  (4.44a) 

^  =  -wh?^t]>  (4-44b) 

$> =**>«)- jrh^*-')'  (4-44c) 

«£?  =  C  (4-44d) 

The  extra  term  in  Eq.(4.41)  now  can  be  written 

<>  (r,  t)  +  3Ap(r,  t)  =  f-a^^  +  3A J  p(r,  t)  =  0,  (4.45) 

where  the  relation  3A(2A  +  vf  =  ua2  in  (4.22)  is  used.  Note  that  this  implies 
A  =  A(n(r,£),T(r,  J)).  The  conditions  for  the  solution  /'"'  are  now  satisfied.  The 
result  (4.36)  can  be  viewed  as  an  extension  of  the  stationary  solution  to  its  local 
form,  just  as  the  local  equilibrium  distribution  is  an  extension  of  the  Maxwellian 
distribution. 
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The  first  order  distribution  function,  /''>,  can  be  obtained  from  Eq.(4.35) 
since  we  now  know  the  expression  for  /'°'.  The  distribution  function,  /(1)  is 
linear  in  the  hydrodynamic  gradients  and  has  the  form 

f^ry,t)  =  xJ^  +  XT,^  +  XuJ-^,  (4.46) 

where  the  Xa<i  are  functions  of  the  velocity  to  be  determined.  Substitution  of 
(4.46)  into  the  (4.35)  provides  a  set  of  coupled  equations  for  these  coefficients, 

-^  -  L(v',a)  +  \6Ui—  +  u\  Xa,k  -  aXUl,k6a,Uy  =  -Ya,k,  (4.47) 

where  a  denotes  n,T,ux,uy,  and  uz,  and  i,k  represents  Cartesian  coordinates. 
The  functions  Yaik  are  given  explicitly  in  terms  of  /'°',  and  their  detailed  form 
can  be  found  in  Appendix  E,  along  with  the  solution  to  (4.47).  Now  the  complete 
normal  solution  near  the  uniform  shear  state  is  constructed  to  order  e.  This  is 
our  primary  result  for  describing  transport  in  states  near  uniform  shear  flow.  It 
has  a  form  similar  to  the  Chapman-Enskog  solution  near  equilibrium.  In  both 
cases  the  deviation  from  the  reference  state  is  proportional  to  the  hydrodynamic 
gradients  and  polynomials  in  the  velocity.  Here,  however,  the  reference  state  is 
not  Gaussian  and  depends  on  the  shear  rate.  In  addition,  the  coefficients,  X0ik 
also  have  a  strong  dependence  on  shear  rates.  Therefore  even  though  this  normal 
solution  is  valid  for  small  spatial  gradients  from  the  reference,  it  can  be  applied 
to  any  finite  shear  rate. 
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4.3     Hydrodynamic  Equations 

The  BGK-Boltzmann  equation  gives  the  hydrodynamic  equations  by  integrat- 
ing the  equations  with  summational  invariants.  They  are  Eqs.(4.8)-(4.10).  These 
are  an  exact  set  of  equations,  but  incomplete.  Once  the  moment  and  heat  fluxes 
are  specified  as  functionals  of  the  hydrodynamic  variables,  a  closed  set  of  equa- 
tions may  be  obtained.  At  the  low  density  near  equilibrium,  they  are  expressed 
by  a  "Newton's  viscosity  law"  for  momentum  flux  and  "Fourier's  heat  law"  for 
heat  flux.  The  hydrodynamic  equations  combined  with  these  fluxes  are  called  the 
Navier-Stokes  equations.  Far  from  equilibrium  state,  no  specific  constitute  equa- 
tions are  known  in  general.  However  by  using  our  normal  solution  from  section 
4.2,  we  can  calculate  these  fluxes  for  the  special  case  of  a  system  near  uniform 
shear  flow.  The  fluxes  can  be  written  with  two  terms, 


U,  =  ^+e^,  (4.48) 


?.  =  <7,(0)  +  «;,(1),  (4.49) 

4r)  =Jdvm(vi  -  Ui)(vi  -  u,)f<-r\  (4.50) 

q«  =  J  rfv^m(v  -  u)2(v  -  u)/<r>.  (4.51) 

t\j  and  q\  represent  the  transport  properties  at  the  local  reference  state  and 
they  are  already  calculated  in  the  previous  section.  The  second  terms,  *[■'  and 
q\     describe  gradients  near  the  reference  state  and  can  be  calculated  from  the 


vhere 
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definition  (4.50)  and  (4.51)  by  using  /'''.  They  are  written  in  terms  of  deviations 
in  the  hydrodynamic  variables  from  the  reference  state, 


and 


*-**;♦<*}•  <«3» 


The  coefficients  of  the  gradients  are  generalized  transport  coefficients  given  by 
ill  =  J  dvm(vi  -  Uj)(vj  -  Uj)Xuuk,  (4.54) 

&j  =  /  dvl-m(v  -  u)2(Vi  -  ut)Xaj,  (4.55) 

where  a  denotes  n  and  T  and  t,  j,  k,  and  /  represent  x,  y,  z.  Their  detailed  calcula- 
tion and  resultant  functions  of  the  shear  rate  can  be  found  in  Appendix  F.  There 
are  many  new  transport  coefficients  which  do  not  exist  for  states  near  equilibrium 
due  in  part  to  the  anisotropy  of  the  reference  state.  Eqs.(4.24)-(4.26)  show  the 
transport  coefficients  at  the  reference  state.  The  heat  flux  was  zero  at  the  state 
because  of  the  definition  of  the  uniform  shear  flow  (uniform  temperature).  Now 
let  us  examine  how  the  thermal  conductivity  changes  near  the  reference  state 
under  the  shear  effect.  The  thermal  conductivity  is  defined  as  the  coefficients 
of  the  temperature  gradients,  so  Q.  •  in  (4.53)  can  be  identified  as  a  generalized 
thermal  conductivity  tensor.  This  is  due  to  the  anisotropy  of  the  reference  state. 
The  diagonal  parts,  Ajj  are  shown  in  Figure  4.2  as  a  function  of  shear  rate. 
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Figure  4.2:  The  diagonal  parts  of  the  generalized  thermal  conductivity,  Xjj(a) 
are  shown  as  a  function  of  a.  The  solid  line  is  for  j  =  y,  the  dashed  line  is  for 
j  =  x,  and  the  dash-dot  line  is  for  j  =  z.  Xjj  is  in  unit  of  kBps/mi/  and  a  is  in 
unit  of  v. 


With  all  the  shear  dependent  generalized  transport  coefficients  known,  the 
new  constitute  equations  are  completed,  and  a  corresponding  the  set  of  closed 
hydrodynamic  equations  is  obtained.  These  equations  are  again  valid  to  second 
order  in  spatial  gradients  relative  to  the  reference  state,  but  there  is  no  restriction 
on  the  value  of  shear  rate  (in  other  words,  there  is  no  restriction  for  how  far  the 
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reference  state  should  be  from  equilibrium  state).    Finally,  the  hydrodynamic 
equations  are 

(Dt  +  Su-  V)n  +  nV-<5u  =  0,  (4.56) 


(A  +  «5u  ■  V)Suk  +  p-1  A[tW  +  tJ*  j£*]  +  „^  =  0,  (4.58) 

where  Dt  =  g  +  Ot^Tj-er,  The  density  equation  is  the  usual  Euler  equation.  At 
the  local  stationary  state  (the  reference  state),  the  temperature  is  constant  due 
to  the  thermostat  (see  Eq.(4.41)).  However,  near  the  uniform  shear  state  the 
temperature  changes  due  to  the  spatial  gradients  of  the  hydrodynamic  variables. 
The  momentum  flux  from  the  reference  state,  <}■  ,  can  be  divided  into  reversible 
and  the  irreversible  parts  as  in  Eq.  (4.44a),  where  the  pressure  gives  the  usual 
Euler  contribution. 


CHAPTER  5 
TWO  TIME  CORRELATION  EQUATION 

The  kinetic  model  equation  for  the  two  time  phase  space  correlation  function 
is  given  by  Eq.(3.38).  In  this  chapter,  this  two  time  correlation  function  is  consid- 
ered for  uniform  shear  flow.  The  equation  for  the  two  time  correlation  function 
has  the  same  structure  as  that  for  the  distribution  function  linearized  around  the 
reference  state.  Hence  the  analysis  of  the  two  time  correlation  function  in  the 
reference  state  provides  a  way  of  studying  both  transport  properties  for  states 
near  uniform  shear  flow  and  the  stability  of  the  reference  state  (a  generalization 
of  Onsager's  assumption  on  the  regression  of  fluctuations  [32,33]).  Similarly,  the 
dynamics  for  fluctuations  of  the  hydrodynamic  fields  in  the  reference  state  is 
the  same  as  the  linearized  hydrodynamic  equations  of  the  previous  chapter  for 
small  perturbations.  To  study  the  transport  properties  and  stability  of  the  states 
near  uniform  shear  flow,  we  solve  the  dispersion  relations  from  these  linearized 
hydrodynamic  equations  and  obtain  the  five  modes  which  govern  the  long  time, 
large  wavelength  dynamics  of  the  system.  For  any  finite  shear  rate  and  suffi- 
ciently long  wavelength,  we  find  that  a  pair  of  propagating  modes  are  unstable 
for  this  flow.  This  new  instability  is  compared  with  a  Monte  Carlo  simulation  of 
the  distribution  function  using  the  Bird  method  [20]  and  shows  excellent  agree- 
ment. To  address  the  question  that  this  instability  might  be  an  anomaly  of  the 
BGK-kinetic  model,  we  study  the  stability  at  small  shear  rates  directly  from  the 
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usual  Navier-Stokes  equations  for  the  hard  sphere  case  at  both  low  and  high  den- 
sities. The  analysis  shows  the  long  wavelength  instability  again.  These  results 
indicate  that  this  instability  is  not  a  consequence  of  the  assumptions  behind  the 
hydrodynamic  calculation  or  use  of  the  low  density  BGK-kinetic  model. 

5.1     Linearized  Hydrodynamic  Equations 

In  Chapter  3,  we  derived  the  kinetic  model  for  the  equation  of  two  time  phase 
correlation  function,  (3.38).  The  formal  solution  (£'  =  0)  to  the  equation  is 
written  in  Fourier  representation, 

C(k,v,t;k',v',0)  =e-'£(k'v'C(k,v,0;k',v',0),  (5.1) 

where  the  operator  C  is  defined  in  Eq.(3.40).  The  time  evolution  operator  in 
the  two  time  correlation  equation  is  the  same  as  that  for  the  BGK-Boltzmann 
equation  linearized  around  the  reference  state,  5f  =  f  —  f„, 

<S/(k,  v,  t)  =  e-tZ{k'v)Sf(k,  v,  0).  (5.2) 

Hence  as  far  as  the  dynamics  is  concerned,  the  two  time  correlation  function  and 
linearized  distribution  function  have  same  time  dependence.  This  is  an  example 
of  a  nonequilibrium  fluctuation-dissipation  relation.  The  two  time  correlation 
function  describes  the  dynamics  of  spontaneous  fluctuations  in  the  reference  state, 
but  it  also  represents  the  linear  response  function  for  small  perturbations  of  the 
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reference  state.  Consequently  by  analyzing  fluctuations  in  the  reference  state,  we 
can  study  how  a  disturbance  of  that  state  decays  or  grows  in  time. 

While  it  is  possible  to  study  directly  the  spectrum  of  C  and  consequently 
all  possible  relaxation  modes,  here  we  focus  instead  on  the  hydrodynamic  spec- 
trum which  is  expected  to  dominate  for  long  times.  The  two  time  matrix  of 
the  hydrodynamic  correlation  function  is  defined  by  velocity  integration  with 
ipa  =  {l,v,t;2}  (see  Eq.(3.10)), 

Ca/3(k,t;k\0)  =  J  dvdV  ija{v)Mv')e~tZ{Kv)C(Kv,0;k',v',0)  (5.3) 

=  J  dvi>a{v)e'tZ(-k'v)  J  dv'^(u')<5(k,v,0;k',v',0),    (5.4) 

and  similarly  the  Fourier  transformation  of  the  linearized  hydrodynamic  variables, 
6ya,  for  deviations  from  the  reference  state  are  given  by 

<5j/a(k,t)  =  |dv^(w)e-'?<k'"><5/(k,v,0).  (5.5) 

Thus  Ca/}(k,  i;k',0)  are  the  response  functions  for  hydrodynamic  perturbations. 
The  advantage  of  this  is  that  we  can  use  the  hydrodynamic  equations  (4.56)- 
(4.58)  which  were  derived  in  the  last  chapter,  to  construct  the  equations  for  both 
Ca0(t)  and  Sya(t). 
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The  deviations  of  the  hydrodynamic  variables  from  uniform  shear  flow  are 
defined  by 

<5n(r,  t)  =  n(r,  t)  -  n„  6T(t,  t)  =  T(r,  t)  -  T„   lu(r, J)  =  u(r,  {)■  (5.6) 

Substitution  of  these  variables  into  the  hydrodynamic  equations  (4.56)-(4.58), 
and  retaining  only  linear  terms  gives 


Dt  Sn  +  nsV  •  <5u  =  0, 


(5.7) 


g2fa         .      d25T 


Dt6uk+p,1 


dn        dr,       [  dT        dr{       1"n-'dridrm 


=  0, 
(5.8) 
+akj8uj  =  0,  (5.9) 


where  s  denotes  a  stationary  value.  These  equations  can  be  written  for  dimension- 
less  hydrodynamic  variables,  5ya,  in  Fourier  representation  form  for  convenience, 


di  +  aijki~dk  )  &Va  +  ^AaB  ~  '  ?  kiBaP<J  ~  S  kjkiDaSji)5yp  =  0,      (5.10) 


where  5ya  is 


Sya  =  <  — , 


3  ST      I   m 


2  t,  '  V  ^ 


-<5u 


(5.11) 
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and  j,  I  denote  x,  y,  z.  The  three  matrices  Aag,  Bag,  Dap  are 


'  0  0  0  0  (T 
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0  0  0  0  0 
000«0 
0  0  0  0  0 
0  0000 


3aft) 
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V' 
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where  j,  /  represents  Cartesian  coordinates.  The  homogeneous  solution  to 


\ 


(5.12) 


(5.13) 


(5.14) 
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equations  (5.10)  can  be  calculated  easily  by  putting  k  =  0, 

5n(t)  =  <5n(0),  (5.15) 

ST(t)  =  ST{0),  (5.16) 

6ux(t)  =  5ux(0)  -  atSuv(0),  (5.17) 

5uv(t)  =  <5u„(0),  (5.18) 

Su,(t)  =  Suz(0).  (5.19) 

These  are  the  expected  results  from  Eqs.(4.37)-(4.39).  Thus  the  homogeneous 
state  is  unstable  to  perturbations  in  5uy,  leading  to  unbounded  linear  change 
in  time.  It  is  expected  that  this  should  be  controlled  at  finite  k  by  exponential 
damping  factors,  e~ak2'  (a  is  a  damping  constant).  This  question  is  addressed  in 
the  next  section. 

5.2     Hydrodynamic  Modes 

For  the  inhomogeneous  case,  to  simplify  the  analysis  we  allow  the  perturba- 
tion only  along  the  velocity  gradient  direction,  i.e.  kx  =  kz  =  0.  In  this  case  the 
linear  hydrodynamic  equations  become  simpler  (i.e.  Dt  — »  <9(),  and  the  hydrody- 
namic modes  can  be  calculated  directly  from  dispersion  relations.  The  linearized 
hydrodynamic  equations  became 


a 

— Sya  +  Fa06yff  =  0,  (5.20) 
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where  Fap  =  Aag  -  ikyB'aB  —  k^D'a0.  The  prime  means  the  case  of  kx  =  kz  =  0. 
Matrices  B't  D'  are  simpler  and  they  are  given  by 


B„a  — 


I 


kBTs 


\ 


0 
0 

w 

5n 


?).VfS(#) 

./ilk 

V3p, 


0  0 


0 


1     0 

2iy      B^yy      0 

0  0  0 
0  0  0 
0        0     0/ 


(5.21) 


DL 


0  0 


\ 


0~l1yx     P~l~lyl 
rs      'xy,s  rs      !yy,s 


n-l-yyy    n~1-^yy 

rs      lxy,s   rs      *yy,s 


0  0 

0  0 

0 
0 


rs     'zy,s  J 


Equation  (5.20)  can  be  solved  by  Laplace  transformation, 


(5.22) 


(zl  +  F)Sy{k,  z)  =  5y{k,t  =  0), 


(5.23) 


vhere 


Then,  the  Sy(k,  z)  is 


5y(k,z)  =  /     dte-"Sy(k,t). 
Jo 


(5.24) 


5ya(k,z)  =  KjSyp(k,t  =  0) 


(5.25) 
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where 

KaB(k,z)  =  [zI  +  F]a0.  (5.26) 

The  determinant  of  the  matrix  K  gives  the  five  simple  hydrodynamic  poles  that 
determine  the  dominant  dynamics  of  the  6ya(k,t)  at  large  t  and  small  k.  At 
equilibrium  (zero  shear  rate),  the  matrix  Kap  becomes  identical  to  the  matrix 
Ba0  in  (3.60)  and  they  give  the  expected  five  hydrodynamic  modes;  two  sound 
modes,  a  heat  mode  and  two  shear  modes  which  were  already  discussed  in  the 
first  example  of  section  3.4.  For  finite  shear  rate,  the  determinant  is  complicated 
and  it  is  very  difficult  to  find  the  dispersion  relations  analytically  .  The  best  way 
to  understand  the  dispersion  relations  for  finite  a  is  to  calculate  each  pole  for  a 
given  value  of  a.  For  convenience,  we  use  dimensionless  variables  for  k,  a  in  units 
of  inverse  mean  free  path,  inverse  mean  free  time  respectively.  Also  since  the  last 
element  in  5ya  is  decoupled  from  other  elements  in  the  matrix  K,  the  dispersion 
relations  are  only  considered  for  four  hydrodynamic  modes.  They  are  obtained 
from  the  equation 

det(Ka0(k,  a,  z))  =  (z  -  zjk, a))(z  -  z2(k,a))(z  -  z3(k,  a))(z  -  z4(k,  a))  =  0. 

(5.27) 
First,  when  the  shear  rate  is  0.1, 


zx  *2$  (0.093  -0.16i)k2/3  +  (-0.75-  1.3 i)  fc4/3  -  0.5 /t2  (5.28) 

z2    =    «J  (5.29) 

z3^-0A9k2,  (5.30) 
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Zi  ^-o.igfc^  +  l.S/c^-O.Sfc2.  (5.31) 

Notice  that  the  real  values  of  z\  and  z-i  are  positive,  which  means  that  the  hy- 
drodynamic  modes  are  growing  in  time  instead  of  decaying  back  to  the  reference 
state.  Also  most  modes  are  not  analytic  functions  of  k,  so  it  is  not  possible  to 
identify  them  simply  in  terms  of  the  equilibrium  hydrodynamic  modes.  The  sys- 
tem is  unstable  for  the  chosen  value  of  the  shear  rate  at  asymptotically  small  k. 
Now  let  us  do  a  numerical  calculation  at  large  k  (k  =  0.2)  with  the  same  value 
of  the  shear  rate  to  see  what  happens  to  the  unstable  modes, 


zi  =  (-0.016  -0.18  i)  (5.32) 

z2  =  z\  (5.33) 

z3  =  -0.02  (5.34) 

z4  =  -0.028  (5.35) 


These  unstable  modes  become  stable  at  this  large  k.  From  this  analysis,  we  may 
speculate  that  there  are  separate  domains  for  stable  and  unstable  dynamics  in  the 
wavevector  and  shear  rate  plane.  The  figure  5.1  shows  the  critical  line  dividing 
the  stable  (above)  and  unstable  (below)  regions  in  the  k  -  a  plane. 

This  instability  is  an  obstacle  in  the  study  of  the  dynamics  of  the  hydrody- 
namic modes.  According  to  the  figure  5.1,  for  any  nonzero  shear  rate  and  small 
k,  there  is  always  an  instability.  It  seems  that  the  dynamics  become  stable  when 


the  wavevector  kvo  >  a,  where  v0  is  a  thermal  velocity,  defined  by  J2kBT/m. 
To  study  the  stable  dynamics  of  the  hydrodynamic  modes,  we  introduce  a  new 
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dimensionless  variable,  x  =  kvo/a,  in  the  matrix  F  of  Eq.(5.20).  This  new  vari- 
able has  a  minimum  number  which  depends  on  the  choice  of  the  shear  rate  value. 
For  example  x  >  0.75  for  a  <  0.3^  to  be  in  the  stable  region.  The  numbers  were 

1.0 


Figure  5.1:  Critical  lines  for  stability  are  shown.  The  BGK-kinetic  model  gives 
the  solid  line  and  the  other  three  lines  are  from  Navier-Stokes  order  for  different 
densities.  Dashed-dot  line  is  for  n*  =  0,  dotted  curve  is  for  n"  =  0.2  and  dashed 
curve  is  for  n"  =  0.4.  All  are  in  units  of  the  mean  free  path  and  mean  free  time. 


determined  from  figure  5.1.  Therefore  the  x  is  not  a  small  variable.  Now  all  the 
shear  rates  in  the  matrix  F  may  be  replaced  by  kvo/x  and  the  matrix  F  becomes 
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a  function  of  k  and  x.  The  Eq.(5.20)  now  can  be  calculated  by  the  perturbation 
method  with  small  k  at  fixed  x  and  the  solution,  5ya(k,  t),  will  show  the  dynamics 
in  the  stable  domain  for  the  hydrodynamic  modes.  The  eigenvalues  of  matrix  F 
are  calculated  from  det(F  -  si)  =  0,  where  s  are  eigenvalues  and  /  is  an  Identity 
matrix.  Also  the  corresponding  eigenvectors,  £,  can  be  obtained 

F(fc,zK<')  =  si(fc,x)C(').  (5.36) 

For  the  convenience  of  the  calculation,  the  dimensionless  variables  are  used  again. 

k*  =  kva/v,    a*  =  a/v,   s*  =  sjv.  (5.37) 

The  eigenvalue  problem  is  solved  for  2-dimensions.  Then  the  dimensionless  eigen- 
values up  to  order  k2  are 


2  ' 

k^_    2k;2 

2    +  5r 

(&.. .  n 

^2      5x 


(5.38) 


+  T-T-  (5-39) 


i\hk'+[--T:^)k'2,  (5.40) 


i\h:k 


'  +  {\-6)k'2-  ^ 


Notice  that  in  this  perturbation  x  is  fixed  at  larger  than  J\.  These  eigenvalues 
are  now  analytic  functions  of  k  in  the  stable  region  and  it  is  possible  to  recognize 
them  as  familiar  sound  modes,  heat  modes  and  a  shear  mode.  The  corresponding 


eigenvectors  are 


c(1)  = 
c(2)  = 


c 


(3) 
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2/3      _j_ 

5'  V5'  ~^/%x 

5'  V  5'      v/53; 

3    /2    /e« 


,0 


52 


,  1 


c 


(4) 


5'  V  5x' 


The  associated  biorthogonal  set  is  denoted  by  {r]^}, 

Q=l 

and  they  are 

rim  =  (— v/10,  0,  -iV&e,  0), 


M/I 


+  — =,  J-,  iv^s,  0 
v/10    ¥5' 


n(3)=|I,/3      »      q.1 


„«) 


2  V  5'  v/Io'    '  2 

2V5'  yio      2 


(5.42) 
(5.43) 
(5.44) 
(5.45) 


(5.46) 


(5.47) 

(5.48) 

(5.49) 
(5.50) 


The  linearized  hydrodynamic  variables,  8ya  can  be  expressed  in  terms  of  these 
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eigenvalues  and  eigenfunctions, 

6ya(k,t)  =  J2e-t3-{k-x)^](k,x)r,{;\k,x)Syg(K0).  (5.51) 

It  is  also  possible  to  give  the  two  time  hydrodynamic  correlation  function,  Eq.(5.4), 
with  these  eigenvalues  and  eigenvectors  because  they  share  the  same  time  evolu- 
tion operator, 

<5o/J(k,t;k',0)  =  ^e-s'<t'I»(C<'»(k,2:)r)«(k,i)Crf(k,0;k',0),  (5.52) 

where  the  CQ/j(k,  0;  k',  0)  is  an  equal  time  correlation  function.  The  discussion 
of  it  will  be  deferred  to  the  next  chapter.  The  two  time  hydrodynamic  correla- 
tion function  can  be  used  to  calculate  the  dynamic  structure  factor,  5(k,tu).  It 
describes  the  density  fluctuation  that  can  be  measured  in  light  scattering  exper- 
iment, 

/oo  „ 

dt  0(k  -  krai„)0(k'  -  kL)C„„(k,  t;  k',  0)e""',  (5.53) 

00 

where  the  0  function  restricts  the  spatial  regions  because  of  the  instabilities.  The 
expression  can  be  written  to  identify  the  modified  Brillouin  and  Rayleigh  peaks, 

S(k,  w)  =  SR(k,  w)  +  5B(k,  w)  +  5s(k,  -w),  (5.54) 
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with 

k2 
SR{k,w)  =  2Re[-iw  +  - — ]-'2Cn(k,  a) + 

2Re[-itu  +  •/-  +  S_l(|  <5n(k,a)  +  ^  <521(k,a))),(5.55) 

SB{k,w)  =  2Re[-i(w  +  c0k)  +  -^ ^-]-1(^<5u(k,a)  +  -^C21(k,  a)),  (5.56) 

pmu     bi/        10  10 

where  f3  =  1/kgT  and  CQ0(k,  a)  are  equal  time  correlation  functions.  Also  c0  is 

a  sound  velocity. 

Now  let  us  examine  more  closely  the  instability  [45]  we  found.     First,  to 

illustrate  how  the  unstable  modes  behave  in  time,  Eq.(5.20)  is  solved  for  given 

numerical  values  of  k  and  a,  which  were  chosen  from  the  unstable  region  in 

figure  5.1.   The  calculations  are  done  with  the  values  of  k"  =  klmIp  =  0.1  and 

a*  =  a/v  =  0.5  and  for  initial  deviations,  5j/(k,  t  =  0)  =  {0,0,  — 0.3j/4o,  »/4o}  are 

chosen  and  y4o  is  constant.    The  detailed  calculation  is  given  in  Appendix  G. 

The  calculation  is  done  for  2-dimensions  because  the  z-component  of  the  velocity 

deviation,  Suz  is  decoupled  from  Eq.(5.20).  This  is  easily  solvable, 

Suz{k,  t)  =  e-^TZ-'SuiQt,  0)  (5.57) 


e  t^r^>J5uz(k,0).  (5.58) 
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Now  the  remaining  hydrodynamic  variables  are 

SUt')  =  2/4o[0.17002e-°'05690m'  -0.00139119e-000320425t" 

(0.83128  cos(0.103845i*)  -  0.111384 sin(0. 103845 i-))e0020532''], 

(5.59) 
%(**)  =  j/4o[1.41137e-°'0569012'"  -0.0730989e~0003204251" 

+(-1.63828 cos(0.103845t*)-  3.68154 sin(0.103845i*))e0020532'"], 

(5.60) 
6y2(t')  =  ij/4o[0.435579e-00569012t"  -0.0422361e-000324251" 

+(-0.393344  cos(0.103845t*)  +  1.09344  sin(0.103845r))e00205321'], 

(5.61) 
SSl{f)  =  M/4o[-0.211397e-°-05690m"  +  0.0307004e-°0032425r 

(0.180697  cos(0.103845  f)  +  0.530312sin(0.103845r))e0020532''], 

(5.62) 

where  t*  is  time  in  units  of  v.  Figure  5.2  shows  the  results  for  5ux(t*)  and 
5uy(t')  as  function  of  time.  It  is  seen  that  there  is  an  unstable  propagating 
mode  of  increasing  amplitude  of  the  maxima  and  minima.  This  is  compared 
with  a  Monte  Carlo  simulation  of  the  BGK-kinetic  equation  using  extended  Bird 
method  [46]  by  Jose  M.  Montanero  and  Andres  Santos.  This  computer  simulation 
is  done  without  any  detailed  hydrodynamic  equation  information.  The  excellent 
agreement  shows  that  the  instability  is  not  a  result  of  the  assumptions  behind  the 
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hydrodynamic  equations.  Even  though  the  computer  simulation  shows  a  good 
agreement,  it  brings  up  the  possibility  that  the  instability  might  be  an  anomaly 
of  the  BGK-kinetic  model,  so  it  might  not  be  a  physical  property.    To  answer 


uy 


Figure  5.2:  Time  evolution  of  (a)5ux(t*)  and  (b)  5uy(t')  for  a*  =  0.5  and  k*  =  0.1. 


this  question,  we  have  calculated  the  dispersion  relations  on  Navier-Stokes  order 
without  kinetic  models.  The  Navier-Stokes  order  means  that  the  gradients  of  the 
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thermodynamic  variables  are  very  small  including  the  shear  rate  as  well,  so  the 
heat  and  momentum  fluxes  are  specified  by  Fourier's  law  and  Newton's  viscosity 
law.  We  consider  this  time  a  fluid  of  hard  spheres,  not  a  Maxwell  molecules 
case,  in  order  to  find  the  effect  of  the  density  dependence  on  the  instability.  The 
equation  of  state  for  hard  spheres  is  calculated  using  Enskog  kinetic  theory  and 
Percus-Yevick  approximation  and  the  transport  coefficients  are  estimated  from 
the  Chapman-Enskog  method  [47—49].  The  linearized  hydrodynamic  equations 
for  ya  are 


+  AN  -  ikBN  +  k2DN    Sya  =  0, 


m 


where  "N"  denotes  Navier-Stokes  order  and  k  =  kv.  The  matrices  are 


(5.63) 


AN 


'  000(V 
0000 
0  0  0a 
000  0 


(5.64) 


on  — 


0 
0 

Po     s  on 

na  dp 
Po  dn 


/2  a_n->   di) 
Zpo±saT 

V  3P0ar 


0 

l\ 

-i/fe*fc 

fi 

0 

0 

0 

% 

(5.65) 


where  e  is  the  internal  energy,  r/  is  the  shear  viscosity  and  po  is  the  pressure  of 
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ideal  gas.  Finally  the  matrix  Dn  is 

/ 


DN 


0      0      0       0 


A  naKB 

0      0      2       0 
p 

0       0       0  \l  +  *l 

3  p         p  J 


(5.66) 


where  k,  p,  are  thermal  conductivity  and  bulk  viscosity,  respectively.   Also  p  = 
n,m  is  the  mass  density. 

The  Explicit  expressions  for  the  equation  of  state  and  transport  coefficients 
are  following.  First,  the  equation  of  the  state  is 


P  =  Po(l+  gflW*), 


(5.67) 


where  po  =  nkgT  and  a  is  the  diameter  of  the  hard  sphere.  The  \  's  the  equi- 
librium pair  correlation  calculated  at  contact  and  is  calculated  from  the  Percus- 
Yevick  approximation, 


x  =  (i- 


■nna3       (imcr3)' 
12 
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)/(!■ 


-)3 


(5.68) 


The  transport  coefficients  are 


14  2 

??  =  -(1  +  —Trna3x)2VB  +  t^ 
X  15  5 


M=  1.002  (-™<t3)2xj7b, 


(5.69) 
(5.70) 
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1  2 

K  =  -(1  +  -nno3x)2KB  +  l*Cv, 
X         5 


(5.71) 


where  fiCv  =  1.002/2.522  (2/37rncr3)2x«;B-   In  these  expressions,  r/B  and  kb  are 
derived  from  the  Boltzmann  equation  for  hard  spheres, 


5    jmkaT     , 

*« = 1016  x  Tsi-r°  ■ 


(5.72) 


ks  =  1-C 


15fcs 
4  m 


(5.73) 


The  dispersion  relations  are  obtained  from  det(FN  -  si)  =0  in  order  to  find 
five  hydrodynamic  modes,  s,(a,  k),  where  Fn  =  Afj  —  ikBfi  +  k2Dx-  They  are 
calculated  for  finite  a  and  asymptotically  small  k  explicitly  to  order  k2. 


I   ->  -<S  +  <  ^^  ~  (J  +  i  ^)2C3*4/3  +  (^  +  C4)^, 


A 


*•     '      v2   ' 

'   2  "            *2 

S2    =    «1, 

S3  -+  Cifc2, 

«4  -»  c2A:2/3  - 

-c3ki/3  +  c<k2, 

s5  ->  V~k2, 

(5.74) 

(5.75) 
(5.76) 
(5.77) 
(5.78) 
(5.79) 


where 


=  kBTnr)dp      2dri. 
2m  p  pdn        dn 


Ci  ■ 


AkBTi^_ 
3m     p 


1/3 


(5.80) 
(5.81) 
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2k rT ,         ,a2      n  dp,  , 
9m  p2      pan 

c,  2k         7n      u  , 

c"  =  T  +  ^—ir  +  r  +  -•  5.83 

3       3n3fcB      3p      p 

These  Ci,C2,C3  and  C4  are  positive  constants  depending  only  on  the  density  and 
temperature  of  the  state.  Since  the  first  two  complex  conjugate  pairs  are  negative 
for  small  k  and  finite  a,  s4  and  s5  are  unstable,  while  the  other  three  modes  are 
stable  for  small  k  and  finite  a.  To  compare  the  results  with  the  ones  from  BGK- 
kinetic  models  (5.28)-(5.31),  the  same  numerical  value  for  shear  rate  (a  =  0.1) 
can  be  used  in  these  s;'s.  They  are 


s_ 

k 


1  ^4  (-0.087  -0.15  i)k2/3  +  (0.80  -1.38  i)ki/3  +  0.6  k2,         (5.84) 
s3  £$  0.20fc2,  (5.85) 

s4  ^4o.174k2/3-1.6k4/3  +  0.6k2.  (5.86) 

Notice  that  since  these  Si  are  from  det(F/v  —  si),  the  negative  means  unstable.  It 
proves  that  the  instability  is  not  a  result  from  the  BGK-kinetic  models.  Figure  5.1 
shows  the  critical  lines  in  the  k-a  plane  separating  stable  (above)  from  unstable 
(below)  dynamics  at  three  different  densities  and  also  shows  the  results  from  the 
BGK-kinetic  model.  When  the  density  is  zero,  the  critical  line  is  matched  with 
the  one  from  the  BGK-kinetic  model  in  the  small  shear  rate  region.  It  shows 
the  same  qualitative  results  for  the  two  different  models.  Since  the  Navier-Stokes 
results  are  calculated  for  small  gradients,  the  curves  for  the  large  shear  rate  region 
are  outside  the  validity  range.  The  "long  wavelength  instability"  is  clearly  seen 
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in  this  figure.  This  instability  gives  a  limited  region  to  study  the  dynamics  of  the 
hydrodynamic  modes.  In  next  section,  the  stable  dynamics  are  discussed  and  the 
two  time  correlation  functions  are  calculated. 


CHAPTER  6 
EQUAL  TIME  CORRELATION  FUNCTION 

In  the  previous  chapter,  we  studied  the  dynamics  of  the  two  time  correla- 
tion functions.  To  complete  the  analysis,  we  need  to  know  the  initial  conditions 
for  those  correlation  functions.  The  initial  condition  is  the  equal  time  correla- 
tion function  and  represents  time  independent  fluctuations  within  the  reference 
state.  While  there  have  been  some  calculations  of  such  fluctuations  in  station- 
ary nonequilibrium  states  [2],  most  applications  were  limited  to  states  near  local 
equilibrium.  Recently,  accurate  experimental  results  [4,29,38]  have  confirmed  the 
predictions  for  these  weakly  nonequilibrium  fluctuations.  Here  we  calculate  equal 
time  correlation  functions  for  a  uniform  shear  flow  state  far  from  equilibrium  us- 
ing the  kinetic  model  equation,  Eq.(3.39),  from  Chapter  3.  The  same  eigenvalue 
method  which  was  developed  in  section  5.2  is  used  and  the  spatial  long  range 
(algebraic  decay)  correlation  is  shown  to  persist  far  from  equilibrium.  However, 
some  limitations  on  the  range  occur  due  to  the  instability  discussed  in  the  last 
section. 

In  Chapter  3,  we  derived  the  kinetic  model  equation  for  the  equal  time  corre- 
lation function,  (3.39).  The  general  solution  to  this  equation  for  the  stationary 
state,  (3.64),  can  be  written  in  Fourier  representation  as 

C,(k1,v1;ka,v2) 
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=  «(ki  +  k2)<5(Vl  -  v2)/s  +  J™  dt  e"(Z'1+^)t7(k1,  vi;  k2,  v2)         (6.1) 
=  <5(k1+k2)<5(v1-v2)/3  +  G,(k1,v1;k2,v2).  (6.2) 

The  second  term,  Gs  is  a  stationary  pair  correlation  function  defined  by  Eq.(3.7). 
We  notice  that  this  pair  correlation  function  has  two  of  the  operators,  C  used 
in  the  two  time  correlation  function,  here  defined  over  functions  of  vl  and  v2. 
Again,  we  limit  attention  to  the  long  wavelength  limit,  for  which  the  equal  time 
matrix  of  the  hydrodynamic  correlations  are  appropriate  (i.e.  Eq(5.3)  with  t  =  0) 

Ga/!(ki,k2)  =  /dv1dv2^a(v1)V/3(v2)C'5(k1,v1;k2,v2),  (6.3) 

=  <5(k1+k2)y"dv^(v)V/,(v)/sW  +  (5^(k1,k2),        (6.4) 

where 

Ga/j(ki,k2)  =  /  dv2dv2)/>Q(vi)i/>/j(v2)G,(ki,v1;k2,v2)  (6.5) 

=  f™dt  /dv1dv2^a(v1)^(v2)e-(Z'+^"7(k1,v1;k2,v2|/s). 

(6.6) 

The  first  term  of  Eq.(6.4)  is  proportional  to  5(ki+k2)  and  represents  short  ranged 
spatial  correlations  as  in  the  equilibrium  case.  The  second  term  of  Eq.(6.4), 
Gap,  represents  the  hydrodynamic  pair  correlation  function.  At  equilibrium, 
the  hydrodynamic  pair  correlation  function,  G0jj  is  zero  because  "/(fe)  in  (6.6) 
vanishes  by  definition  in  (3.36). 
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However,  at  nonequilibrium,  7  is  no  longer  zero  and  provides  new  nonequi- 
librium  phenomena.  At  the  states  near  equilibrium,  the  density  fluctuations  for 
uniform  shear  flow  already  have  been  calculated  in  the  second  example  of  section 
3.4.  They  showed  spatial  algebraic  decay  of  the  fluctuations  in  (3.76).  For  states 
far  from  equilibrium  (large  shear  rate),  the  hydrodynamic  equal  time  correlation 
function  also  can  be  analyzed  by  using  the  eigenvalues  and  eigenvectors  of  £j  and 
£2.  However,  there  is  a  limitation  for  the  spatial  range  because  of  the  instability 
discussed  in  the  last  chapter.  Then  the  hydrodynamic  eigenvalues  should  be  cal- 
culated for  limited  region  k  >  kc(a)  (see  Figure  5.1).  Since  the  operators  £j  and 
£2  are  the  same  as  those  used  for  the  two  time  correlation  functions  in  the  last 
chapter,  the  same  hydrodynamic  eigenvalues,  s,  in  (5.38)-(5.41)  and  the  same 
eigenvectors,  C(l'  in  (5.42)-(5.45)  can  be  used  for  the  analysis  of  the  pair  corre- 
lation function.  Then  the  hydrodynamic  pair  correlation  (6.6)  can  be  written  in 
terms  of  these  eigenvalues  and  eigenfunctions  and  the  time  integral  performed  to 
give 

Go^k,, k2)  =  £  ]  Cii)(k1,x)Cy)(k2, »)fl*(k,kfc x),   (6.7) 

where  x  =  kv0/a  is  the  variable  used  in  the  section  5.2.  The  matrix  //IJ(k1,k2,a;) 
is  defined  by 

fl*(k,,ks,*) 

=  Ef?7)(ki,i)^')(k2^)/dvl(iv2^(v1)^(v2)7(k1,v1;k2,v2)  (6.8) 


86 

=  «S(ki  +  k2)  £ „W(k!, x)rj«'(k2, 1)1/  / dv^(v)^(v)(/.(v)  -  &(▼)). 

10  J 

(6.9) 

The  second  line  of  the  equation  is  obtained  using  ■yBGK  in  (3.36)  and  the  elements 
of  H'i  are  given  in  Appendix  H. 

Let  us  now  specialize  Eq.(6.7)  to  the  density  fluctuations  for  uniform  shear 
flow,  i.e.  a  =  0  =  1, 

_Ve   2(t»  +  ,V5)  23]  +  _A_3 

where  fc  and  a  are  used  in  units  of  v/vo  and  1/  for  simplicity.  This  expression 
for  the  density  fluctuation  is  restricted  to  k  >  v/fa,  and  according  to  Figure 
5.1  this  k  and  a  relation  is  valid  for  a  <  0.5.  When  a  =  0  (equilibrium  case), 
Gn  becomes  zero  which  is  the  expected  equilibrium  result,  and  for  finite  a  the 
correlation  function  shows  long  ranged  spatial  correlations. 


CHAPTER  7 
CONCLUSIONS 

The  objective  in  this  thesis  has  been  to  study  transport  and  fluctuations 
for  states  far  from  equilibrium  for  the  special  nonequilibrium  states  of  uniform 
shear  flow.  Dynamical  properties  of  states  far  from  equilibrium  are  not  well 
understood  due  to  their  complexity  and  technical  difficulties  with  the  formal 
theories  of  nonequilibrium  statistical  mechanics.  However  at  low  density,  kinetic 
theory  provides  a  controlled  formulation  of  this  problem.  It  consists  of  three 
kinetic  equations:  the  Boltzmann  equation,  and  related  equations  for  the  two 
time  and  equal  time  correlations.  They  describe  all  states  from  equilibrium  to 
far  from  equilibrium.  There  are  still  difficulties  for  practical  applications,  except 
for  cases  at  or  close  to  equilibrium. 

In  this  thesis,  we  constructed  kinetic  models  corresponding  to  these  kinetic 
equations  in  order  to  resolve  practical  difficulties  for  system  far  from  equilibrium. 
The  simple  BGK  model  for  the  Boltzmann  equation  has  been  extended  in  a  self- 
consistent  way  to  apply  to  kinetic  equations  for  fluctuations  as  well.  All  of  the 
important  structural  relationships  of  the  three  kinetic  equations  are  preserved  by 
this  model.  If  the  solution  to  nonlinear  BGK-Boltzmann  equation  is  obtained, 
the  linear  equations  for  the  two  and  equal  time  correlations  can  be  solved.  These 
kinetic  model  equations  were  applied  to  uniform  shear  flow  to  calculate  transport 
properties  and  fluctuations  for  large  shear  rate.  First,  the  BGK-Boltzmann 
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equation  was  solved  for  uniform  shear  flow  and  the  normal  solution  was  obtained 
near  uniform  shear  flow  by  using  a  variant  of  the  Chapman-Enskog  approximation 
method.  By  "variant"  we  mean  that  in  this  expansion  for  uniform  shear  flow,  the 
localized  form  of  the  stationary  solution  is  used  as  a  reference  function  instead 
of  the  local  equilibrium  distribution  function.  This  normal  solution  consists  of 
hydrodynamic  variables  and  their  spatial  gradients, 

/  =  /(°)(n(r,  t),T{t,  t),  u(r,  i))  +  /(1)(Vn,  VT,  V<5u).  (7.1) 

This  function  then  was  used  to  calculate  the  irreversible  momentum  and  heat 
fluxes  in  terms  of  the  hydrodynamic  gradients, 

where  "*"  denotes  the  irreversible  part  of  these  fluxes.  The  coefficients  of  these 
gradients  are  transport  coefficients.  Since  these  coefficients  are  calculated  near 
the  stationary  state  rather  than  at  the  local  equilibrium  state,  there  are  many 
new  coefficients  (f]l,  Qj,  (Li),  which  characterize  the  transport  properties  of 
the  system  far  from  equilibrium  (i.e.,  for  large  shear  rate).  With  these  new 
expressions  for  the  fluxes,  a  closed  set  of  the  generalized  hydrodynamic  equations 
was  derived  near  uniform  shear  flow. 

The  equation  of  the  two  time  correlation  function  for  hydrodynamic  fluctu- 
ations in  uniform  shear  flow  has  the  same  dynamics  as  that  for  the  linearized 


hydrodynamic  equations  around  the  reference  state.  We  studied  the  excitations 
(hydrodynamic  modes)  of  the  uniform  shear  flow  by  linearizing  these  hydrody- 
namic equations.  It  was  verified  that  there  exist  hydrodynamic  modes  even  far 
from  equilibrium  (i.e.  modes  that  vanish  as  k  — >  0).  However,  we  also  found  that 
there  is  an  instability  for  long  wavelengths  and  any  finite  shear  rate.  This  in- 
stability is  compared  with  computer  simulation  results  (by  J.  Montanero  and  A. 
Santos)  and  showed  very  good  agreement  between  theory  and  simulation  during 
the  initial  stages  (see  Figure  5.2).  Because  of  the  linear  analysis  in  theory,  the 
large  time  behavior  could  not  be  studied,  but  the  simulation  has  been  performed 
for  longer  times  and  shows  continuing  growth  of  the  hydrodynamic  variables. 
For  the  final  stage  of  this  instability,  the  simulation  predicts  a  homogeneous 
state  with  a  higher  temperature.  It  is  likely  that  the  long  wavelength  instability 
has  not  been  seen  in  previous  molecular  dynamics  simulations  due  to  the  finite 
system  sizes  considered. 

Finally  the  equal  time  correlations  were  investigated.  The  long  ranged  corre- 
lation in  the  density  fluctuations  was  calculated  for  small  and  large  shear  rate. 
For  states  near  equilibrium  (small  shear  rate),  the  spatial  algebraic  decay  of  the 
fluctuations  has  been  known  for  some  time.  Similar  results  were  obtained  here  for 
large  shear  rate  as  well.  However  there  is  a  spatial  restriction  this  time  because 
of  the  long  wave  length  instability.  The  effect  of  this  instability  on  both  equal 
and  two  time  correlations  requires  further  study. 

In  this  thesis,  we  applied  the  special  nonequilibrium  state  of  uniform  shear 
flow  to  "model"  kinetic  equations  in  order  to  study  transport  and  fluctuations. 


These  model  equations  also  can  be  used  for  steady  heat  flow  or  combined  heat 
and  momentum  fluxes,  for  which  exact  solutions  to  the  BGK-Boltzmann  equa- 
tion are  known  [50,51].  This  will  be  deferred  to  future  work.  Also,  it  remains 
to  understand  the  nature  of  the  instability  and  its  consequences.  The  directions 
of  the  perturbations  can  be  generalized  (e.g.,  kx  #  0)  for  analyzing  the  hydrody- 
namic  modes.  Then,  we  can  obtain  a  more  complete  description  of  long  range 
correlations  (e.g.,  we  can  have  a  directionally  dependent  long  range  correlation). 
By  using  a  BGK-Enskog  equation,  the  system  can  be  extended  to  higher  densities 
where  we  can  study  the  instability  for  conditions  more  suitable  to  molecular  dy- 
namics simulation.  At  high  density,  a  large  k  hydrodynamic  instability  has  been 
proposed  [52, 53]  and  shown  to  occur  at  shear  rates  similar  to  those  observed  in 
the  molecular  dynamics  simulations  [9].  Further  study  is  required  to  relate  the 
mechanism  in  the  large  k  hydrodynamic  equations  to  that  described  in  this  thesis 
at  small  k. 


APPENDIX  A 
DENSITY  FLUCTUATIONS 

In  Chapter  2,  density  fluctuations  are  studied  and  long  range  spatial  correla- 
tions are  verified  for  small  shear  rate.  In  this  appendix,  a  few  steps  to  reach  the 
results  will  be  shown.  The  nonequilibrium  part  of  density  fluctuations  (2.34)  is 
written  by  using  the  analysis  from  Dufty  [27]  and  Groome,  Gubbins  and  Dufty  [28] 
as 

lim  (6n{r)6n(r')Mxy{t))e  = 

t—*oo 

-7^3  /  dkl7^eiMr_r)  /     dT(P*<  nJ.[S(k,v),3(-k,v)]) 
(27r)J  J         oky  Jo 

(A.1) 
where  (X,  Y)  is  defined  by 

(X,Y)=fdpf0XY,  (A.2) 

and 

a(k,v)  =  e-l(n,-'"-v)l.  (A.3) 

The  a  function  depends  on  the  inhomogeneous  linearized  Boltzmann  operator 
nl  —  ik  ■  v  and  can  be  represented  in  terms  of  its  proper  eigenvalues  and 
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eigenvectors.    To  solve  equation  (A.l),  first  write  the  function  a  in  terms  of 
combination  of  eigenvectors  of  nl  —  ik  •  v 


(nl  -  ik  •  v)ga  =  saga, 


(A.4) 


where  sa,  ga  are  eigenvalues,  eigenvectors  respectively.  This  equation  is  solved 
by  perturbation  theory  for  small  wavevectors.  The  dominant  contribution  for 
long  range  correlations  arises  from  the  hydrodynamic  spectrum  (eigenvalues  that 
vanishes  as  k  — >  0).  The  hydrodynamic  eigenvalues  and  eigenvectors  can  be  found 
in  the  literature  [54],  so  only  the  results  are  given: 


12  f    e         ik.f    „      1   ,,       lA.f      , 


93 

54 


^§fh+>'k-r^-^kjf-^ 

m    f .  ik f _.  \ 

I  m    /  ik- 


kg!        £       p 


(A.5) 
(A.6) 
(A.7) 
(A.8) 


with  eigenvalues 


s1|2  =  ±i  k  c  +  r&2,  53  =  s4  =  r]/(nm)k2,  s5  =  DTk2 


(A.9) 


Here  the  wavevector  is  k  =  kk  and  the  orthogonal  directions  to  k  are  given  by 
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the  unit  vectors  {ei,e2}.  Also  e  =  |m»2  and  S,  T  are  solutions  to 


IS,  =  5„  /7Ij  =  %  (A.10) 

where  5j  and  T;_,  are 

Si  =  ("imw2  -  -k8TJ  Vt,   Ttj  =  m  (viVj  -  -8,jV2^  .  (A.ll) 

The  solutions,  S  and  T  exist  due  to  the  orthogonality  of  Si  and  T^  to  the  sum- 
mational invariants.  In  Eqs.(A.5)  and  (A. 8),  8  =  2Km/5kBrj  is  a  dimensionless 
parameter  and  k,  jj  are  the  thermal  conductivity  and  shear  viscosity,  respectively. 
In  the  eigenvalues,  c  =  J5kBT/Zm  is  the  sound  velocity,  and  T  and  Dt  are  the 
sound  damping  constant  and  thermal  diffusivity.  The  normalization  of  eigenvec- 
tors is  chosen  such  that  (ha,  gp)  =  5ap,  where  ha(=  g*)  are  complex  conjugates  of 
ga.  Now  the  function  5  is  expressed  in  terms  of  these  hydrodynamic  eigenvectors, 

afovj-^ce— <*>rSo.  (A.12) 

a=l 

Consequently  Js[a(k),  a(-k)}  can  be  written 

Js[5(k,v),5(-k,v)]  = 
Clc2e-2r*v(7s[5l,  hi]  +  Js[g2,  h2})  +  4e-1B^rJ,\gt,  h5] 

+c1cse-^k<+n2+D^(J3[gu  hs]  +  J3[g5,  h2\) 
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+c2c5e'^-ikc+rk2+D^2\j5[92,  h,]  +  JJ*,  fcj),  (A.13) 


where 


3/2       ifc   /    m      f 1 


C5  =  -</|.  (A.15) 


Use  of  this  expression  for  /,[5(fc),  5(— k)]  in  Eq.(A.l)  shows  that  the  first  term  on 
the  right  side  of  Eq.(A.l)  vanishes  because  the  term  (pI,nJs[5(k,  v),  5(-k,v)]) 
becomes  zero  due  to  the  momentum  conservation.  Then  Eq.(A.l)  is 

\im{6n{r)6n(r')Mxy{t))e 

(A.16) 

r  roo 

(A.17) 
Substitution  of  Eq.(A.13)  in  the  integral  of  the  right  side  of  (A.17)  leads  to 


(%y,nJs[a(k),a{-k)})  =  -2clc2(3n{Txy,Txy)kxkylk'1  e~2Tk  T 

1.2/ „2„-2Tiifcc    ,    „2„2Tiic\„-2Trfc2 


+0n(Txy,Txy)kxky/k2(c2e-2T'kc  +  c22e2T'k')e- 
-2ik-^==[2{Txy,{t  +  3/2T)Txy)kxky/k2  -  (Txy,eTJ}) 
+m(Txy,  {vxSy  +  vySx))kxky/k2  -  m(Txy,  (vySx  +  vxSy))kxky/k2 
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2li,„2^-2nfcc        „2„2Tikc\„-2Trk' 


-2*7/2(1  -  9/2)(Txy,Txy)kxky/k2](cle-2^  -  c22e2^)e 

o3/2 

+2i  k-r — \m{Txy,  {vxSy  +  vySx))kxky/k2 
5m 

-m^,,,  (^5„  +  vySx))kxky/k2  -  2ev(Txy,Txy)kxky/k2 

-2(Txy,  (e  +  3/2T)Txy)kxky/k2  +  (T^eTjj)  +  5(Txy,Txy)kxky/k2} 

x(c1c5e-"A:':  -  c2c5eT'kc)e-T^+D^k2 .  (A.18) 

After  performing  time  integral  and  inverting  the  Fourier  transformation,  the  lead- 
ing term  of  the  nonequilibrium  part  of  the  density  fluctuation  for  large  \r  -r'\  is 
found  to  be 


}un(5n{r6n(r>)Mxy(t))e  =  l^^faQ  +  0[r~2}.  (A.19) 


Therefore  the  density  fluctuation  for  uniform  shear  flow  is 


(Sn(rSn(r'))ne  ~  {8n(v)5n(T'))e  -  aA  JL^_^_^} .  (A.20) 


APPENDIX  B 
DERIVATION  OF  KINETIC  EQUATIONS 

In  this  appendix,  the  three  kinetic  equations  (3.11)-(3.13)  are  derived  from 
an  approximation  of  BBGKY  hierarchy  by  an  expansion  in  powers  of  a  low  den- 
sity parameter.  The  expansion  parameter  is  obtained  by  scaling  space  and  time 
with  characteristic  scales.  Since  we  are  interested  in  the  phenomena  that  change 
appreciably  over  scales  such  as  mean  free  path  (E0  S  l/no2)  and  mean  free  time 
(i0  =  io/vg),  they  are  used  as  characteristic  units  to  form  a  dimensionless  hierar- 
chy. Here,  n  is  the  density,  a  is  the  force  range,  v0  =  (2fcsT/Tn)'/2  is  the  average 
thermal  velocity,  and  kg,  T  and  m  are  Boltzmann  constant,  temperature  and 
the  particle  mass.  With  these  scalings,  the  reduced  distribution  function  /(s)  is 
scaled  with  n' /vl". 

To  simplify  the  analysis,  we  consider  a  gas  of  hard  spheres,  but  the  same 
results  can  be  obtained  for  general  short  ranged  forces  by  more  elaborate  methods 
[19,30].  Since  the  hard  sphere  interaction  is  singular,  the  forces  are  replaced  by 
binary  scattering  operators  [31].  The  resulting  dimensionless  BBGKY  hierarchy 
has  the  form 


s  +  5>«-V,-a»£T(i,i)   /«(*.—,»„«) 

01        1=1  Kj  J 

=  £/(%,+,  r(M+i)/<»+l>{j/i,.  •-,},,+!,*).  (b.i) 
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T(i,j)  is  an  operator  describing  the  scattering  of  the  pair  of  particles  {i,j}, 

T(i,j)  =  /*  6(r  •  glJ)f  ■  fe  (<5(r,;  -  ar)bt]  -  *(r«  +  erf)) ,  (B.2) 

where  df  denotes  a  two  dimensional  solid  angle  integration  over  the  sphere  for 
particle  at  contact.  Also  gy  =  Vj  -  Vj  is  the  relative  velocity,  f  is  a  unit  vector 
along  the  ri3 ,  and  by  is  the  scattering  operator  defined  for  any  function  X(Y{,  v,) 

by 

b,rY(v„vy)=X(v;,v;.).  (B.3) 

The  scattered  velocities  are  vj  =  vt  -  (gy  ■  f)f  and  v'j  =  v3  +  (gy  ■  f )f. 

Equations  (B.l)  and  (B.2)  depend  on  the  single  dimensionless  parameter  a  = 
a/£0  =  nc3,  the  ratio  of  the  force  range  to  the  mean  free  path.  This  parameter  is 
small  at  low  density  and  can  be  used  as  an  expansion  parameter  for  the  solutions 
to  Eq.(B.l).  The  s-particle  reduced  distribution  is  expanded  in  powers  of  a, 

/W(yi,  ■  ■  ■ ,  V.,  t)  =  £  a2"/M(2/i,  -,*,*).  (B.4) 

n=0 

Substituting  Eq.(B.4)  into  the  Eq.  (B.l)  and  equating  coefficients  of  a2"  gives  the 
kinetic  equations  determining  the  coefficients  /^s'.  The  dependence  of  a  in  the 
operator  T(i,  j)  is  neglected  because  the  distances  of  interest  are  large  compared 
to  the  force  range,  and  the  position  difference  between  colliding  partners  can  be 
neglected.  Consequently  T(i,j)  simplifies  to 

T(i,  j)  =  6(rt])  J  di  6(r  •  glJ)f  ■  to(^  -  1).  (B.5) 
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To  zeroth  order  in  a  the  solutions  are  found  to  be 


/n*.  •••.  v.,  *)«n /«%.*).  (B-6) 

where  /q     is  the  solution  to  the  kinetic  equation 

(j^  +  v!  -V,j  /PW)  =  / 'dlteT(l,2)/i1)(y1,0^1,(lft,*)-  (B-7) 

The  binary  collision  operator,  the  right  side  of  the  Eq.(B.7),  is  recognized  as  the 
nonlinear  Boltzmann  operator  for  hard  spheres,  so  /0(     is  indeed  the  solution  to 
the  Boltzmann  equation. 
The  solution  to  order  a2  is 

/i(2)(yi,J/2,*)=/(S1)(2/i,t)/1(1)(y2,t)  +  ^1)(!/2,i)/l(1)(2/1,t)  +  G1(y1,y2;t).     (B.8) 

The  s-particle  distribution  function  /{  can  be  determined  in  terms  of  single 
particle  function  /g  ,  /}  and  pair  function  G\,  but  will  not  be  considered  here. 
The  functions  /}     and  G\  are  the  solutions  to 


d 

^  +  v1-V1 


Ai)  f[l\yut)  =  J  dy2  T(l,  2)G1(y1,  y2,  t),  (B.9) 


f^+v1-V1-A1+v2.V2-A2JG1(!/l,!/2,t)  =  r(l,2)/(S1)(!/1,i)^1)(y2,*)., 

(B.10) 
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Here  Ai  is  the  operator  defined  over  functions  of  yt  by 

A,X(y,)  =  J  dy2T(\,2)Ul>l)(yut)X(y2)  +  X(yi)f^(y2,t)).  (B.ll) 

This  a  expansion  gives  the  exact  solution  to  the  entire  hierarchy  up  to  order  a2. 
The  reduced  distribution  functions  are  determined  as  a  sum  of  products  of  single 
particle  functions  /^(j/i,*),  f\  (j/i><)  and  the  pair  function  Gi(yi,y2,t).  From 
the  definition  of  single  time  pair  correlation  G(yi,y2,t)  in  (3.7)  and  the  zeroth 
order  solution  (B.6),  it  is  easily  shown  that  the  pair  function  Gi(j/i,lft,i)  is  a 
leading  contribution  to  the  density  expansion  of  G(y\,y2,  t).  Therefore  in  the  low 
density  limit,  G(yi,y2,t)  is  the  solution  to  Eq.(B.lO). 

The  analysis  of  correlations  at  two  times  can  be  accomplished  in  a  similar 
way.  Define  the  two  time  distribution  function 

h^l\yu---,ya,t;y\t')  =  (i^(yi,---,ys,t)}^{y',t')).  (B.12) 

The  two  time  correlation  function  in  (3.6)  is  related  to  /i'2>  by 

C(y,t :  y',t')  =  h™  (y ,  t;  y' ,  t')  -  fm(y,t)f^(y',t').  (B.13) 

The  two  time  distribution  function  /i<s'  obeys  the  same  BBGKY  hierarchy  (B.l) 
for  t  >  t'.  When  a  power  expansion  in  a  is  applied,  the  solution  up  to  order  a2  is 

4S+1V-  •  •  ,ys,t;y',t')  =  jfty.O/J'W  ■  ■  .«**).  (B.14) 
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M2)(yi,t;y',t')  =  /^)(y1,i)/i(1)(y'.i')  +  /r)fei,t)^I>(y',*'),  (B.15) 

and  for  s  >  2 

h[s+l\yu-  ■  ■  ,ys,t;y',t)  =  f{0s)(yu-  ■  ■  ,ys,t)A1](y',t') 
+As)(yu-  ■  ■  ,ys,t)f^(y',f)  +  i:tlfl)l){yJ,t)Cl(yl,f)y',t'), 

(B.16) 

where  C\  obeys  the  equation 

f^+v-V-AJc1(3/,«;y',i')  =  0>  (B.17) 

for  i  >  J'.  Substitution  of  (B.16)  into  (B.13)  shows  that  Ci(y,t;y',t')  is  the  first 
nonvanishing  contribution  to  C(y,  t;  y',  £'),  so  that  C(y,  t;  y',  t')  satisfies  Eq.(B.17) 
at  the  low  density  limit.  This  completes  the  derivation  of  all  low  density  kinetic 
equations  for  transport  and  pair  correlations  at  one  and  two  times. 


APPENDIX  C 
DERIVATION  OF  KINETIC  MODEL 

In  this  appendix,  the  operator  ABGK  (3.31)  is  derived  in  detail.  The  model  for 
the  linear  operator  A  is  obtained  from  the  functional  relation  to  the  Boltzmann 
collision  operator, 

abgk%)  =  JdVl  (-JL-j**W{t)))  %l} 

=  -v(h(y)-Jdy1h(y1)jj^-T)ft(y,t)\.  (C.l) 

Here  JBGK  =  —v(f  -  fi)  is  used.  The  functional  dependence  of  fi  on  /  occurs 
only  through  the  the  functions  za.  These  functions  are  defined  by 

J dvf^Y)f(y,t)  =  J dvfa(y)Uy,t),  (C.2) 

where  ipa  =  {l,v,ii2}.  This  definition  assumes  that  the  BGK  model  retains  the 
exact  conservation  laws.  Therefore  the  functional  derivative  form  of  fi  on  /  is 
written  as 

-Ji{y,t)  =  fi(y,t)Mv)r(,,  ,,za(r,t).  (c.3) 


5f(yuty™   >       "«■    >™    >Sf(yi,t) 
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Take  the  functional  derivative  of  each  side  of  (C.2)  to  get 

5 


S(t  -  ii)ipavi  =  J  dvipa(v)  — — xft(V,t) 


S 


J  dvipa{v)fe(y,  t)i>0(v)———zg(T,  t) 


where  gap  is  defined  by 

ga0(r,  t)  =  Jdv  ft(y,  t)7/-Q(v)^(v).  (C.5) 

Inversion  of  (C.4)  for  the  functional  derivatives  gives 

^-^Q(r,  t)  =  8(v  -  r,)  (£}(rf  t)M^)-  (C.6) 

Put  this  form  into  (C.3)  to  get 

47(^lj /«(»■')  =  ft(y,t)S(r-ri)M^9^t)M^i)-  (C7) 

This  result  is  used  with  Eq.(C.l)  to  give  the  BGK  model  for  operator  A 

ABGK%)  =  -1/(1  -  P)h{y),  (C.8) 

where  P(t)  is  defined  in  Eq.(3.32),  and  is  a  projection  operator  along  the 
summational  invariants,  ipa. 


APPENDIX  D 
EQUILIBRIUM  EXCITATION 

The  dynamics  of  correlations  at  an  equilibrium  state  is  determined  from  the 
matrix  Bap  =  5ap  —  Iap  ,  Eq.  (3.60),  which  determines  both  the  hydrodynamic 
and  kinetic  poles.  The  matrix  Iag  is  obtained  from  the  Fourier-Laplace  transform 
of  (3.51),  specialized  to  the  equilibrium  state 

Iaff(k,z)  =  I/^1|dv/,(v)fl(k,Z|v)^(v)^(v),  (D.l) 

where  R(k,z\v)  is  given  by  (3.59).  Further  analysis  is  simplified  by  choosing  the 
collisional  invariants  i>a(v)  to  be  an  orthonormal  set,  </>a(v),  constructed  from 
linear  combinations  of  those  in  (3.22), 


,   mv         3\       /  m  I  m  I  m 

u  ^3  law5 "  2 J '  foUl'  fe'2'  fe"3 1 


where  V\  =  k  •  v  is  the  component  of  v  along  k,  and  V2  and  V3  are  the  orthogonal 
transverse  components.  With  this  choice,  gap  =  n5ap  and  (D.l)  reduces  to 

laS(k,z)  =  ~Aag{x).  (D.3) 

kv0 

Here  the  dimensionless  variable  x  =  (z  +  i/)/ikv0  is  used,  and  v0  =  (2kBT/m)1'2 
is  the  thermal  velocity.  Aa/}  is  a  symmetric  matrix  and  its  elements  are  given  in 
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terms  of  the  complex  probability  integral, 


Z(x) 


/wJ-t 


tfr 


dv 


exp  —v 


by 


(Aa0) 


Z       {\y*\x  +  (a?-\)Z\     V2(l  +  xZ)0  0 


.4121[Z  +  (|)^(^- 

-  iMd 

\Z2xAi2 

0  0 

A«                            A23 

V^x^is 

0  0 

0                      0 

0 

Z  0 

0                      0 

0 

0  Z 

(D.4) 


(D.5) 


Since  Ia0  is  analytic  for  Re(z)  j&  —v,  the  singularities  of  Ba0(k,  z)  leading  to 
hydrodynamic  modes  are  associated  with  values  of  z  for  which 


det  (Ba0)  =  det  (6a0  -  Ia0)  =  0, 


(D.6) 


and  vanish  as  k  — >  0.  This  equation  is  solved  numerically  to  look  for  five  hydro- 
dynamic  dispersion  relations,  z  =  z(k),  such  that  z(k  =  0)  =  0.  In  fact  due  to  the 
structure  of  the  matrix  Aap,  Eq.(D.6)  factorizes  into  two  independent  equations. 


det(Ba/3)  .  |(1 -  l-^-)2(M(k,  z)  +  iN(k,  z))  =  0.  (D.7) 


The  first  equation  gives  a  twofold  real  degenerate  solution,  zs(k),  and  there  is 
a  maximum  value  of  k,  k*  =  y/Trv/v0,  where  the  first  equation  no  longer  has 
a  solution.    When  k  =  fcj,  the  value  of  z(k")  =  -v  .    This  branch  represents 
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the  shear  modes  defined  for  k  <  k*.  To  lowest  order  in  k,  one  gets  the  familiar 
Navier-Stokes  expression  zs  =  -r/fc2,  where  r\  =  p/v  is  the  kinematic  viscosity. 

The  other  equation,  M  +  iN  =  0,  leads  to  two  solutions  which  are  complex 
conjugates  and  to  another  real  solution.  These  are  the  two  sounds  mode,  z±,  and 
the  heat  mode,  Zh ,  respectively.  To  Navier-Stokes  order,  z±  =  ±ick  -  Vk2  and 
zh  =  -Dxk2,  where  c  is  the  sound  velocity,  T  is  the  sound  damping  constant,  and 
Dt  is  the  thermal  diffusivity.  For  the  BGK  model,  (j  =  T  =  DT,  in  agreement 
with  what  is  observed  in  Fig. (3.1).  The  sound  modes  and  the  heat  mode  also 
present  limiting  values  of  k.  From  our  numerical  analysis  we  have  found 
k±  ~  1.85i//«o  and  k*H  ~  1.91z//uo,  corresponding  to  z±  =  zh  =  —v. 


APPENDIX  E 
FIRST  ORDER  APPROXIMATION 

In  this  appendix,  the  first  order  equation,  (4.35)  is  solved.  Eq.(4.35)  is 

(AW  +  V  •  V)  /<°)  =  -  (^  ~  W> a>  +  XSu>iy.  +  u)  f®>  (E-1) 

where  Dt  =  Tr+<*y*V^r-  The  terms  on  the  left  side  of  Eq.(E.l)  can  be  written 
as 

(flj»>  +  v- .  V)/(°)  =  -£  ^(A(1)  +  V  •  V)ya,  (E.2) 

where  ya  denotes  {n(i,  £),T(r,  J),  <5u(r,  £)}.  Those  time  derivatives  can  be  re- 
placed by  first  order  spatial  gradients  of  hydrodynamic  variables  by  using  the 
corresponding  hydrodynamic  equations  obtained  from  Eqs.(4.8),  (4.9)  and  (4.10) 
to  first  order  in  the  gradients 


£>,(1)n  +  V  •  rathi  =  0,  (E.3) 


\nkB(D^  +  «5u  ■  V)T  +  atM  +  4°'^  +  V  ■  q(0)  =  0,  (E.4) 

i  '  OTj 

Dit1)6ul+5u-V6ul  +  p-l-^-  =  0.  (E.5) 
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Use  of  these  equations  into  Eq.(E.2)  gives 


(^+v.v)/^yB,|+^+^.  (E.6) 


The  coefficients,  y0|j,  are  given  by 

where  i,  /,  fc  denote  s,  j/,  z.  The  derivative  of  /'°'  is  simply  obtained  by 

■+(g)y>-u         (E.10) 

CE.11) 

(E.12) 
Here  the  functions  ha  are  defined  with  time  integration, 

hn  =  J     drre-^+^VK^A.^-r)^-^)),  (E.13) 

hT  =  Jg    dre-^^e^lA^-r)^  -  5u3)}2  ft(e*  .^(-^(v)  -  fa,)),     (E.14) 


9/'°'  _  /(»)   ,   (dv\  1, 
9n 


5T 

3 

~~2T 

/»  + 

m 
2£BT2 

a/to)  _ 

m 

ffcs,. 
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Kt  =  ^    dr  e-"r+3^e2"rgk(V,  r)/<(eATA,3(-r)(U;  -  8uj)),  (E.15) 

where  ^(v'.t)  in  the  hUk  are 

Si  =  (»i  -  SuT)  +  at(v'y  -  Suy), 

gy  =  at(v'x  -  6ux)  +  (1  +  a2t2){v'y  -  6uy), 

gz=v'z-6uz.  (E.16) 

Now  the  left  side  of  Eq.(E.l)  can  be  written  in  terms  of  known  coefficients,  Yaj 
and  hydrodynamic  gradients.  The  right  side  of  Eq.(E.l)  can  be  written  in  terms 
of  hydrodynamic  gradients  and  unknown  coefficients,  Xaj,  which  was  defined  in 
(4.46).  By  equating  the  coefficients  of  the  hydrodynamic  gradients  in  both  side 
of  equation  finally  gives  the  equations  for  Xa<k. 

-£■  -  L{v\  a)  +  X6u,—  +  v\  Xa,k  -  aXUxtk5a^  =  -Ya,k,  (E.17) 

where  a  =  n,T,ux,uy,uz.  The  last  term  on  the  left  side  of  Eq.(E.17)  exists  due 
to  the  time  dependence  of  the  homogeneous  solution  of  5ux(t)  in  Eq.(5.17).  This 
equation  has  a  form  similar  to  Eq.  (4.34)  and  the  general  solution  can  be  written 
as 

Xa,k  -  jf     dre-"Te3ArHtt,UvX0l,t(eATAu(-r)(W;.  -  5u3)) 

-Y^i^Aiji-r^-SuM-  (E.18) 


APPENDIX  F 
TRANSPORT  COEFFICIENTS 

In  this  appendix,  the  transport  coefficients  defined  in  Eqs.(4.54)  and  (4.55) 
are  calculated  by  using  the  solution  to  Eq.(4.47).  First,  the  coefficients  from 
momentum  flux  are  considered, 

7im  =  -  /  dv'mCiCjXu,tm,  (F.l) 


where  c  =  v'  -  i5u.  As  an  example,  the  yf*  is  analyzed.  Since  the  7^  is  defined 


with  XUxiV,  the  equation  for  XUliV  is  considered  to  solve 


a<°)  d 

-£■  -  L(v\  a)  +  ^ut—  +  v  )  Xu„y  =  -Yu„y,  (F.2) 


where  the  operator  L  is  defined  in  (4.17).  Substitution  of  the  solution,  XUx:y 
(E.18)  into  the  definition  of  f%>,  (F.l)  gives 

7*  =  -  JdV mcxcy  J°°  dt  e-^e^Y^y^Aiji-ty,).  (F.3) 

The  t|*  can  be  calculated  with  couple  of  change  of  variables  on  velocity  variables. 
First  v'  — >  v'  +  <Su  and  then  v[  — >  e~Xiv[,  and  finally  do  v't  — Y  Aij(t)v'j.  Then  7^ 
is 


=  -  /     die  ("+2A)'m  f  dv  (vx  -  atvy)vy 
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3nfcB   &T  {tly+aimlly}+  d5uxy 


«(«8?+»o( 


2i/q  lOi^a 


(2A  +  y)3      3(2A  +  y)3/       (2A  +  i/)2 


(F.4) 


jP.  (F.5) 


Use  of  the  result  of  (ffi  in  (4.44b)  and  the  relation  of  A  and  shear  rate  a  (4.22) 
gives  the  final  value  of  the  7JJ, 


'IV  /n^     ,    ..\').(c\     1    .A"* 


(F.I 


,I!/      (2A  +  y)2(6A  +  i/) 
The  same  method  is  used  for  all  the  transport  coefficients  and  they  are  as  follows: 
2av(v2  -  18i/A-72A2) 


lilb  = 


t2/p 


3  (v  +  2  A)3  (v  +  6  A)2 
(y  +  2A)2(i/  +  6A) 


(F.7) 


(F.I 


TS/p  =  0 


(F.9) 


7J»/p  =  ->6  +  28  y5  A  +  417;/1  A2  +  3087  1/3  A3  +  11622  v2  A4 

+21420 1/ A5  +  15336  A6)/((i/  +  2  A)3  (1/  +  3  A)2  (1/  +  6  A)2) 


(F.10) 


iSlp : 


a  1/  (5 v*  +  76  ^3  A  +  408  v2  A2  +  912  v  A3  +  720  A4) 
3(y  +  2A)5(i/  +  3A)(y  +  6A) 


(F.ll) 


Ill 


iS/P  =  °  (F.12) 


iS/P  =  "ft/P  =  0  (F.13) 


2a  A  (6A  -  5i>) 
^  -  (2A  +  ,V(6A  +  ^  (F'14) 


7« /P  =  (5  ^  +  119  ^7  A  +  1107 1/6  A2  +  2361  vh  A3  -  26388  vA  A4 

-192096 p3  A5  -  512784  u2  A6  -  597456  v  A7  -  233280  A8)/ 

(3^(i/  +  2A)3(t/  +  3A)3(^  +  6A)2)  (F.15) 


„      _  4  a  (2 1/5  +  26  i/4  A  +  166  v3  A2  +  549  v1  A3  +  852 1/  A"  +  468  A5) 
7l»/p_  3(i/  +  2A)4(j/  +  3A)2(i/  +  6A)  (        > 


iTJP  =  0  (F.17) 

■y**/p  =  (a(33i/10  +  1084i/9A  +  15968i/8A2  +  115260i/7A3 

+345843i/6A4  -  431136i/5A5  -  6653448i/4A6  -  21439728i/3A67 

-32523120i/2A8  -  22243912i/A9  -  4199040A10))  / 

(9^(2A  +  ^5(3A  +  ^)4(6A  +  y)2)  (F.18) 

7^/p  =  (2j/6  +  15;/ A  -  16i/A2  -  135i/3A3  +  222f2A4 

+252i/A5  -  1368A6)  /(3i/(2A  +  i/)4(3A  +  i/)2)  (F.19) 
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1?JV  =  0  (F.20) 

iff/P  =  %IV  =  0  (F.21) 

7f*/p  =  (2(^7  +  87i/A  +  1521i/5A2  +  10497;/ A3  +  34344z/4A4 
+57528//A5 

+56592i/A6  +  37584A7))/ 
(3i/(2A  +  i/)3(3A  +  ^)2(6A  +  i/)2)  (F.22) 


2(,3  +  24,2A  +  72,A2) 
7ll/P_     3(*/  +  2A)>  +  6A)2  l        j 


.  2 

7™/p  =  3(J/  +  2A)3(i/  +  6A)  (F'24) 


iX/P  =  0  (F-25) 

7™/p  =  (2  a  i/  (2  v6  +  48 1/5  A  +  278 i/4  A2  +  3  y3  A3 
-4422  v2  A4  -  129961/ A5  -  11448  A6))/ 
(3(i/  +  2A)!>  +  3A)2(i/  +  6A)2)  (F.26) 


_  2  u  (-2  i/3  -  11 1/2  A  -  6 1/  A2  +  24  A3) 
7w,/P"  3(i/  +  2A)4(i>  +  3A)  (R27) 
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7»/p  =  0  (F.28) 

7g/P  =  TJJ/P  =  °  (F-29) 


s        _  2(i/5  +  12i/4A  -  36f3A2  -  582i/2A3  -  1512z/A"  -  1080A5) 
7"'P  ~~  3(2A  +  iv)3(3A  +  j/)(6A  +  ^2  (        ' 


2^  +  24A,  +  72A2) 
7**/P~    3(2A  +  i/)*(6A  +  i/)*  l        J 


7^/p~     3(2A  +  i>)3(6A  +  i>)  (F'32) 


TJI/P  =  0  (F.33) 

7jJ/p  =  (2a(-!/7  -  5j/6A  -  4I1/A2  -  573i/4A3  -  2982i/3A4 
-5508i/2A5  -  648i/A6  +  5184A7))/ 
(3(2A  +  z/)5(3A  +  t>)2(6A  +  i/)2)  (F.34) 


2;/(6A  +  ^)(^2  +  14A^  +  30A2) 
3(2A  +  ^)4(3A  +  y) 


-Tyy/P  =      V    Z^lZrVr::,        '  (F-35) 


7y"/P  =  0  (F.36) 


114 


l"Jp  =  l"jp  =  0  (F.37) 


_  4t/(-i/4  -  bav\  -  U\v*  +  6aA2  -  72AV  -  168^A3  -  144A4) 
7sz/P  "  3(2A  +  y)4(6A  +  y)2  (F-38) 


iTJp  =  Jiy/P  =  0  (F.39) 


*fr— Hx+31  (F-40) 


7*i7p  =  79i*/P  =  7*?/p  =  0  (F.41) 

Tj'/p  =  -(j/6  +  7Ai/5  +  48f4A2  +  216^3A3  +  447j/2A4 

+399i/A5  +  144A6)/(iv(A  +  i/)(2A  +  i/)2(3A  +  i/3))  (F.42) 


„.        2va{v2  +  5Ai/  +  7A2)  „     . 

^/P=     (2A  +  ,)4(3A  +  ,)  (R43) 


iZlv  =  0  (F.44) 

7£/p  =  1%IP  =  lllIP  =  0  (F.45) 

7„Sx/P  =  7^/p  =  0  (F.46) 


115 


"y'/P  -  -Jv-^y  P*0 


a(z/2  +  5Xu  +  8A2)  ,„      . 


2A 


T»;/p  =  0  (F.50) 

-J '—&.,  =  (-(5  va  +  186  v1  A  +  2215 i/  A2  +  13480 1/5  A3  +  49144 1/4  A4 
+112992 1?  A5  +  159552  i?  A6  +  122688  v  A7  +  36864  A8))/ 
(2i/(iy  +  2A)2(i/  +  3A)3(i/  +  4A)3)  (F.51) 


i,bj 


-pi=^  =  (a  (9 1/9  +  108 1/8  A  -  641 1/7  A2  -  16802  i>6  A3  -  116528  v5  A4 
-419184 1/4  A5 -879808 1/3  A6 
-1096896 1/2  A7  -  772992  v  A8  -  248832  A9))/ 
(2i/(i/  +  2A)3(j/  +  3A)3(y  +  4A)4)  (F.52) 


?_«    =  -(5^13+al-11A  +  278I/12A  +  26aI/10A2  +  7865i/uA2 


nklT^T-x 


+289  a  y9  A3  +  1 1 1924  j/10  A3  +  1780  a ;/  A4  +  926254 1/9  A4 
+6560  a  i/  A5  +  4918220 1/8  A5  +  14464  oi/6A6  +  17797024 1/7  A6 
+17664  avs\7  +  45657536  i/6  A7  +  9216  a  i/4  A8  +  85373952  i/s  A8 
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+  118065024  j/4  A9  +  119017728 1/3  A10  +  81077760 i/2  A11 

+30633984  v  A12  +  3538944  A13)/ 

(2 1/2  (u  +  2  A)3  (i/  +  3  A)4  (i/  +  4  A)5)  (F.53) 

-S=S „  =  (14 a v10  +  317 a  v9  A  +  vw  A  +  3078  a  i/8  A2  +  24  j/9  A2  +  17835  av1  \ 
ikiT     * 


+245  v8  A3  +  75678  a  i>6  A4  +  1378  v1  A4  +  268688  a  v5  A5  +  4608  i/6  A5 


nk2BT 


+776528  a  i/  A6  +  9152  i/5  A6  +  1594080  av3  \7  +  9984  p4  A7 
+2017152  a  i/2  A8  +  4608  y3  A8  +  1323648  a  f  A9  +  304128  a  A10)/ 
(2i/(y  +  2A)3(y  +  3A)4(y  +  4A)4)  (F.54) 

f*^  =  (-5y7  -  42i/A  +  203i/5A2  +  3566</4A3  +  16364^3A4 
+35376i/2A5  +  37536t/ A6  +  16128A7)/ 
(2i/(2A  +  y)(3A  +  i/)3(4A  +  y)3)  (F.55) 


m 


,  =y£*  .  =  (a  (-r>13  -  401  i>12  A  -  10705  u11  A2  -  135383  z/10  A3  -  1053900  y9  A 


-5568818  is8  A5  -  20669200  v1  A6  -  53491128 1/6  A7  -  91320528 i/5  A8 
-85889792  z/  A9  -  4286976  v3  A10  +  86142080  v2  A11  +  86512896  P  A12 
+29505024A13  +  24883;0A14))/ 
((»/  +  2  A)4  (*  +  3  A)5  (//  +  4  A)5  (k  +  6  A))  (F.56) 


I  n^L  ~  (A  (29  "    +  739  "    A  +  9717  v  A  +  78521 "  A 
kBT 

+383576  v1  A4  +  994366 i/6  A5  +  404636  ub  A6  -  5601664  u"  A7 
-18060048 1/3  A8  -  26083104  u2  A9  -  18698112  v  A10  -  5336064  A11))/ 


fc|T2 
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[v(v  +  2  A)2  (i/  +  3  A)5  (1/  +  4  A)4  (1/  +  6  A))  (F.57) 

g,  =  (2A(4j/11  +  113i/0A  +  1218^9A2  +  6285i/8A3  +  13216;/ A4 
-17468i>6A5  -  170096j/5A6  -  460452i/4A7  -  848880i/3A8 
-1461168i/2A9  -  1912896fA10  -  1133568A'l))/ 
(y(2A  +  i/)2(3A  +  i/)5(4A  +  k)4(6A  +  u))  (F.58) 


APPENDIX  G 
THE  HYDRODYNAMIC  MODES 

In  this  appendix,  the  calculation  of  the  hydrodynamic  modes,  5ya(t),  for  spe- 
cial case  a*  =  0.5  and  k"  =  0.1  is  done  from  Eq.(5.25)  to  compare  the  behavior 
of  the  instability  with  the  results  from  simulations.  The  simulations  were  per- 
formed using  a  thermostat  with  constant  A,  rather  than  A  =  A(n(r,i),T(r,  £))  as 
considered  here.  To  account  for  this,  the  linearized  equation  for  the  temperature 
in  (5.8)  has  two  additional  terms  that  couple  the  temperature  to  itself  and  the 
density.  These  arise  since  there  is  no  complete  cancelation  of  viscous  heating 
to  lowest  order  in  the  spatial  gradients.  For  convenience,  we  use  dimensionless 
variables;  i.e., 

- k/u,  (G.l) 

m 

A*  =  X/u,    a'  =  ajv.  (G.2) 

From  now  on,  we  considered  only  two  dimension  case  (x,y  dimension),  since  z- 
component  is  not  coupled  to  other  components  so  it  gives  a  result  (5.57).  Also 
we  drop  the  *  sign  from  now  on.  Then  the  determinant  of  the  matrix  K  becomes 
a  fourth  order  polynomial  function  for  z, 

det[Ka0(k  =  0.1,  z)}  =  2.04302  10"6  +  0.000666015  z  +  0.00891949  z2 

+0.0190414  z3  +  z4,  (G.3) 
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and  the  solutions  for  det(K)  =  0  are 

Zl  =  -0.0569012,    z2  =  -0.00320425,   z3  =  0.020532  -  JO.  103845,   z4  =  z'3. 

(G.4) 
The  dynamics  of  6y~a(i)  in  (5.20)  is  determined  from  these  poles.  The  first  two 
poles  are  negative,  so  their  modes  are  stable,  whereas  the  complex  conjugate  pair 
z3  and  24  represents  unstable  modes  since  the  real  parts  of  the  poles  are  positive. 
It  is  easy  to  see  this  situation  if  we  write  down  the  expression  of  5ya(t)  with  these 
poles.  The  Eq.(5.25)  can  be  written 

r»«<*  =  °-^  =  (2  -  gg^^p  -  h)  SUk  =  aM  =  0)-  (a5) 

where  K~p  =  det?L ,  and  the  inverse  Laplace  transformation  of  5ya  gives  complete 
expression  of  5ya(t)  in  terms  of  initial  conditions, 

It  is  clearly  seen  that  23  and  24  are  responsible  for  instability.  The  matrix  elements 
ofNap  inEq.(G.6)  are 

JVn  =  0.000936892  +  0.00366708  z  +  0.0190414  z2  +  z3,  (G.7a) 

Nn  =  -0.0000114814-0.00361146  2,  (G.7b) 

N13  =  0.000131178  i  -  0.00011728  i  z,  (G.7c) 

Nu  =  0.0000607019  i  +  0.0010161  i  z  +  0.0707107  i  z2,  (G.7d) 
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N21  =  -0.00129203  -  0.00625125  z  -  0.0962409  z2,  (G.7e) 

N22  =  0.0000140618  +  0.00443992  z  +  0.00711708  z2  +  z\  (G.7f) 

N23  =  -0.000181528  i  -  0.000287276 1  z  -  0.0367099  i  z2,  (G.7g) 

AT24  =  -0.0000435735  i  +  0.01160132  z  +  0.0698092  zz2,  (G.7h) 

N3l  =  0.00209147 i-  0.0348971iz-  0.0103748 iz2,  (G.7i) 

N32  =  -0.000081186! -0.0254766JZ- 0.0225905 iz2,  (G7j) 

JV33  =  -0.000284938  +  0.00887352  z  +  0.0165958  z2  +  z3,  (G.7k) 

JV34  =  -0.000144986  -  0.00361283  z  -  0.496755  z2,  (G.71) 

Af41  =  0.0000288927 i  -  0.00383077  i  z  +  0.0742802  i  z2,  (G.7m) 

N42  =  0.000162372  i  z  +  0.0510737  i  z2,  (G.7n) 

AT43  =  ,  0.00185513  z  -  0.00165859  z2,  (G.7o) 

yV44  =  0.000858455  z  +  0.0143699  z2  +  z3.  (G.7p) 

Now  we  need  to  determine  the  initial  conditions  for  5ya(t)  in  order  to  see  the 
dynamics  of  the  hydrodynamic  modes  with  instability.  We  chose  the  initial  con- 
dition such  that  6uy(t)  is  dominated  mostly  by  unstable  modes.  The  results  are 
shown  in  section  5.2  and  the  figure  5.2  shows  the  results  for  Sux(t)  and  5uy(t)  as  a 
function  of  time  and  they  are  in  good  agreement  with  the  results  from  simulations. 


APPENDIX  H 
CORRELATION  FUNCTIONS 

The  matrix  H''  is  defined  in  (6.8).  The  elements  of  H'*  =  W*  and  are  given 


by 


1A 


Hn  = 

5x2     "     , 
2X  +  u 

Hl2  = 

-Hn, 

Hu  = 

2V5*(2A  +  ^' 

Hu  = 

-H13, 

H22  = 

jW, 

H23  = 

\[h~H" 

H2i  = 

{&  +  H" 

H33  = 

'o        x 

10         2(2A  +  i/) 

rrU  _ 

Ln^         A 

if44  =  #3 


(H.l) 
(H.2) 
(H.3) 
(H.4) 
(H.5) 

(H.6) 

(H.7) 

(H.8) 

(H.9) 
(H.10) 


where  O  =  i  \ P&.  +    28"«2,  4.    n™' 

wnere  V  -  6  [    4A+„  +  (4a+„)3  t  (4a+./)5 
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